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ABSTRACT 


A  computer  code  was  developed  for  predicting  the  convective  ond 
radiative  heat  transfer  to  the  surface  of  a  flat  plate  with  hypersonic 
boundary  layer  flow  and  with  injection  of  absorbing,  emitting  ond 
scattering  particles  due  to  ablation.  The  program  was  used  to  study  the 
heat  transfer  rate  as  a  function  of  wall  and  freestream  temperatures  and 
Mach  Number  of  the  flow.  The  results  were  also  obtained  to  study  the 
effect  of  injection  of  absorbing  and  scattering  particles  in  the  boundary 
layer  on  the  heat  flux  as  a  function  of  wall  ond  freestream  temperatures, 
Mach  Numbers  and  particle  radiative  properties.  It  is  seen  that,  if  there 
is  strong  external  incident  radiation,  the  injection  of  scattering  particles 
reduces  the  radiative  flux  by  as  much  as  75%  depending  on  the  scattering 
coefficient  of  the  particles  and  the  particle  density  in  the  flow. 

A  review  of  the  published  literature  shows  that  the  data  of  the 
radiative  properties  of  particles  is  locking.  Also  the  data  on  the 
radiative  properties  of  modern  composite  material  surfaces  is  not 


available. 


NOMENCLATURE 


A  tabulated  coefficient  dependent  on  the  particle  size 
parameter  and  refractive  index 

Defined  in  equation  (39) 

Effective  absorptivity  in  the  presence  of  scattering 
Defined  in  equation  (40) 

Diffuse  surface  radiosity 

Mass  fraction  of  injected  species 

Coefficient  of  Diffusion  of  Specie  j  into  specie  i 

Diameter  of  spherical  particles 

Black  body  emissive  power 

Hemispherical  radiant  flux 

Nth  order  exponential  integral  function 


Defined  in  equation  (10) 
Stagnation  enthalpy 


Radiation  intensity 

Change  in  intensity  due  to  absorption,  emission,  or  scattering 


Coefficient  defined  by  K  *  o(a+28) 

Prandtl's  mixing  length 
Coefficient  defined  in  equation  (26) 

Real  part  of  the  complex  refractive  index  m(l-ix) 


Pressure 


Prandtl  number 


Nth  order  Legendre  polynomial 
Convective  heat  flux  at  the  plate  surface 


Radiative  heat  flux 


Total  heat  flux  due  to  convection  and  radiation 
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R 

Volumetric  density  of  spontaneous  radiation 

Rf 

Reflectivity 

Re 

X 

Local  Reynolds  number 

s 

Scattering  function 

T 

Temperature 

u 

x-  component  of  velocity 

V 

y-  component  of  velocity 

X 

Distance  parallel  to  the  plate  surface 

y 

Distance  normal  to  the  plate  surface 

% 

>. 

i 

z 

Particle  size  parameter 

.> 

Greek 

Letters 

Q 

Absorption  coefficient 

> 

. 

6 

Extinction  coefficient 

. 

y 

Scattering  coefficient 

6 

Coefficient  defined  in  equation  (28) 

« 

•> 

6 

X 

Boundary  layer  thickness 

e 

Polar  angle  between  the  incident  ray  and  the  scattered  ray 

M 

Cosine  of  the  angle  between  the  direction  of  propagation  and 
the  x-  axis 

“eff 

Effective  viscosity 

P 

Density 

f 

Azimuthal  angle 

«' 

w 

Solid  angle 

O'O' 

Tabulated  parameters  relating  and  8^  of  the  particles 

l 

T 

Optical  coordinate 

1 

» 
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T 

o 

Optical  thickness 
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in 

Subscripts 


a  Pertaining  to  absorption 

e  Pertaining  to  emission 

eff  Effective  property 

s  Pertaining  to  scattering  of  beam  initially  in  the  direction 

of  propagation 

s'  Pertaining  to  scattering  of  beams  from  all  directions  into 

the  direction  of  propagation 

x  Pertaining  to  location  x 

w  Pertaining  to  wall 

•  Pertaining  to  freestream 

A  Monochromatic  or  wavelength  dependent  property 


INTRODUCTION 


Prediction  of  radiative  and  convective  heat  transfer  to  the  surface  of 
a  hypersonic  space  vehicle  is  very  important  in  the  design  and  development 
of  such  vehicles.  The  predicted  heat  flux  at  the  surface  gives  an 
Important  boundary  condition  for  the  analysis  of  heat  transfer  through  the 
skin  to  the  inside  of  the  vehicle.  However,  this  is  a  complicated  problem 
especially  when  the  medium  of  the  boundary  layer  participates  in  the 
radiative  heat  transfer.  If  the  aerodynamic  heating  is  so  high  that  the 
vehicle  surface  material  ablates,  the  problem  becomes  even  more 
complicated,  since  the  mass  injection  due  to  ablated  material  not  only 
affects  the  flow  in  the  boundary  layer,  it  also  affects  the  radiative 
transfer  because  the  particles  of  injected  mass  absorb,  emit  and  scatter 
thermal  radiation.  The  ablation  phenomenon  can  be  used  os  a  tool  for 
structural  cooling.  And  if  Intense  incident  radiation  is  expected,  such  os 
from  a  nuclear  blast  or  a  laser  weapon,  the  injected  particles  can  be  used 
to  scatter  and  reflect  the  incident  radiation  away  from  the  vehicle 
surface.  This  can  also  be  achieved  by  injecting  particles  with  high  back 
scattering  coefficients,  into  the  boundary  layer. 

This  study  was  undertaken  to  develop  a  computer  code  for  predicting 
convective  and  radiative  heat  transfer  to  the  surface  with  hypersonic 
boundary  layer  flow  and  with  injection  of  absorbing  emitting  and  scattering 
particles  into  the  boundary  layer.  This  was  done  by  incorporating  a 
radiative  heat  transfer  analysis  developed  by  the  author  [1]  into  a 
turbulent  boundary  layer  code  developed  by  Patankar  and  Spalding  [2] .  The 
mass  injection  into  the  boundary  layer  was  added  into  the  computer  code  to 
study  the  effect  of  ablation  [3]. 

The  computer  program  developed  in  this  project  was  used  to  predict  the 
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radiative  and  convective  heat  transfer  rotes  to  a  flat  surface  from  the 
boundary  layer  and  from  an  external  radiative  source.  The  heot  transfer 
rates  were  studied  as  a  function  of  wall  and  free  stream  temperatures  and 
Mach  Number  of  the  flow.  Then  the  results  were  also  obtained  to  study  the 
effect  of  injection  of  absorbing  and  scattering  portlcles  in  the  boundary 
layer  on  this  heot  flux  to  the  surface,  as  a  function  of  wall  and  free 

stream  temperatures,  Mach  Numbers  and  particle  radiative  properties.  It  is 

2 

seen  that  if  there  is  strong  incident  radiation  (we  assumed  1,000  Btu/ft  - 
sec.)  the  injection  of  scattering  particles  reduces  the  radiative  heat 
transfer  by  as  much  as  1*  to  755*  depending  on  the  scattering  coefficient  of 
the  particles  and  the  particle  density  (or  volume  fraction).  This  happens 
without  much  change  in  the  temperatures  in  the  boundary  layer  because  the 
reduction  in  the  incident  radiation  is  mostly  due  to  reflection  to  the 
outside  rather  than  absorption.  If  there  is  no  incident  radiation  from  an 
external  source  and  the  thermal  radiation  is  mainly  due  to  high  free-stream 
temperatures,  the  injection  of  scattering  portlcles  is  not  effective  in 
reducing  the  radiative  heat  flux,  because  the  scattering  particles 
themselves  emit  thermal  radiation  at  the  boundary  layer  temperatures  and 
offset  any  reduction  by  scattering.  However,  the  injection  from  the  sur¬ 
face  reduces  the  temperature  gradient  at  the  surface  thereby  reducing  the 
convective  heat  transfer  and  hence  the  total  heat  transfer  to  the  surface. 

The  results  from  this  study  give  a  boundary  condition  for  heat  flux  to 
the  surface  of  a  vehicle  which  can  be  used  to  predict  the  heat  transfer  and 
temperatures  inside  a  composite  material  wall  of  a  space  vehicle.  The 
results  for  conditions  other  than  those  presented,  can  be  obtained  by 
running  the  computer  code  developed  in  this  study.  It  must,  however,  be 
pointed  out  that  this  study  has  some  limitations  which  must  be  Improved  and 
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it  raises  some  questions  that  must  be  addressed.  A  limitation  of  this 
development  Is  the  way  It  models  ablation.  The  ablation  process  Is  modeled 
simply  as  Injection  of  mass  Into  the  boundary  layer,  without  ony 
consideration  of  the  enthalpy  of  ablation  process  and  the  various  species 
obtained  and  the  surface  modification  due  to  ablation.  A  review  of  the 
literature  conducted  in  this  study  shows  that  the  data  on  the  radiative 
properties  (absorption  and  scattering  coefficients)  of  particles  Is  almost 
non-existent.  Also,  the  data  on  the  radiative  properties  of  modern 
composite  material  surfaces  Is  not  available.  Another  question  that  thl6 
study  raises  is  how  the  surface  ablation  of  the  composite  material  affects 
the  surface  roughness  which  then  affects  the  drag  and  the  convective  heat 
transfer  to  the  surface.  It  16  hoped  that  the  deficiencies  In  the  computer 
modeling  and  the  other  question  raised  In  this  study  can  be  tackled  In  the 


Phase  II  of  the  project. 


THEORETICAL  DEVELOPMENT 


In  order  to  analyze  the  radiative  and  convective  heat  tronsfer  in  a 
turbulent  boundary  layer  over  a  wall  with  surface  ablation,  the  problem  is 
formulated  as  follows: 

Figure  1  shows  the  flow  situation  and  the  choice  of  coordinate  system. 
The  problem  is  formulated  as  a  steady  turbulent  boundary  layer  flow  of  air 
past  a  flat  plate  ot  zero  angle  of  attack.  Effect  of  ablation  is  modeled 
os  mass  Injection  from  the  wall  into  the  boundary  layer.  There  is  assumed 
to  be  no  shock  layer  present  in  order  to  keep  it  purely  boundary  layer 
problem.  However,  output  from  a  shock  analysis  could  be  used  as  input  to 
this  boundary  layer  problem  in  order  to  effect  a  solution  to  the  entire 
flow  field. 

The  governing  equations  for  steady,  two-dimensional,  radiating, 
turbulent  flow  are: 

(a)  the  conservation  of  mass 

57  U)  ♦  ^  (p*)  -  0  <’> 

(b)  the  conservation  of  X-momentum 
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ay  'weff .  ay'  dx 


(c)  the  conservotion  of  energy 


(2) 
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(3) 


the  concentration  equation  for  the  Jth  specie: 


pu 


(4) 


h 


where  u.  v.  p,  H  and  C  represent  the  time  averaged  values  of  the  fluctuate 

r* 

lng  turbulent  quantities.  M  and  Pr  are  called  “effective"  transport 

eff  eff 

properties  and  refer  to  the  turbulent  and  lominar  contributions  to  shear 
and  heat  flux  respectively.  Expressions  for  these  quantities  are  given  in 
Appendix  A.  Since  the  flow  is  steady  ot  zero  angle  of  attack,  it  is  quite 
easy  to  show  that 


—  =  0  for  the  flat  plate 
dx 

The  lost  term  in  the  equation  (3),  divQ  ,  is  called  the  radiation 

R 

source  term.  In  a  high  speed,  viscous  boundary  layer,  locally  the  tem¬ 
perature  gradients  in  the  Y-direction  will  be  for  greater  than  the  tem¬ 
perature  gradients  in  the  X-directlon.  For  this  reason,  a  one  dimensional 
radiative  transport  regime  is  assumed. 

Therefore, 


divqR 


3^  %  1 


(5) 


This  term  can  be  calculated  using  a  simplified  analysis  developed  by 
Goswami  and  Vachon  [l]  and  described  in  the  next  section. 
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RADIATIVE  HEAT  TRANSFER  ANALYSIS 


This  section  describes  the  generol  equotlons  of  rodiotive  heot  tronsfer 
in  absorbing,  emitting  and  scattering  medium  and  olso  describes  the 
development  of  a  simplified  analysis  of  the  above  problem.  The  simplified 
analysis  was  first  developed  by  Goswaml  and  Vachon  [1]  and  was  called  as 
effective  absorptivity  analysis.  Since  then,  the  approach  has  become 
accepted  in  the  scientific  community  and  is  now  known  as  "SCALING"  [4-6]. 

A.  General  Equations  of  Radiative  Transfer 


In  order  to  formulate  the  equation  for  the  rodiotive  flux  vector  in  an 
absorbing,  emitting  and  scattering  medium,  it  is  necessary  to  develop  the 
expression  for  the  Intensity  of  radiation,  1^ ,  within  the  medium.  To  this 
end,  the  development  of  Sparrow  and  Cess  [7]  is  referred  to. 

Referring  to  Figure  2,  the  intensity  of  radiation  within  the  medium, 

■f 

lx(x,0),  can  be  considered  to  be  composed  of  two  contributions,  i.e.,  Ix(x,6) 
directed  in  the  positive  x  direction  and  I^(x,9)  directed  in  the  negative  x 


direction,  such  that: 


lA(x.e)  «  I*(x,e)  -  l"(x,e) 


To  find  expressions  for  1^  and  1^,  a  radiation  balance  is  performed  on 

a  volume  element  of  the  medium  as  shown  in  Figure  2.  As  a  monochromatic 

♦ 

beam  of  intensity  I^(x,  €&  passes  through  the  volume  element,  its 

+ 

intensity  is  changed  by  an  amount  dl  .  This  change  consists  of : 

+  * 

1.  Attenuation  of  lx  due  to  absorption  within  the  element  by  an 
amount  J^.a; 

♦ 

2.  Attenuation  of  1^  due  to  scattering  by  suspended  porticles  in 
the  medium  or  by  the  molecules  of  the  medium  by  on  amount  J^.s; 

+ 

3.  Augmentation  of  1^  by  on  amount  J^.s!  os  0  reSult  of  energy 
scattered  into  the  beam  from  other  directions; 

♦ 

4.  Augmentation  of  Ix  due  to  emission  from  the  elementol  volume 
by  an  amount  Jx , e . 


Vr  < 


V  *  • v,1 


FIGURE  2.  RADIATION  BALANCE  IN  AN  ELEMENTAL  VOLUME 


Therefore,  dl^  can  be  written  as: 


dI=J  -  J  ♦  J  .  ♦  J 

A  A, a  A,S  A  ,S  A,e 


(7) 


Substituting  the  expressions  for  the  term  on  the  right  hand  side  of 

+ 

equation  (7),  Sparrow  and  Cess  [  7]  obtained  the  equation  for  1^  as 


<11  (t  ,p,4>)  .  a 

-3- - Ij  bi 
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(8) 


In  the  same  manner,  they  obtained  the  equation  for  1^  as: 


,u,*)  _  a  Y, 

w  - - -  «=  -I  (T,g,4)  +  — -  eK,(T)  +  7T7-  G,(tJ 


dt. 
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where 


2n  1 


f_y  dp'd#' 


(10) 


^A  ln  thG  above  e<1uation  is  called  the  extinction  coefficient 
(B^=a^+Y^)  and  is  called  the  dimensionless  optical  thickness  given  by 


T 


A 


(11) 


Equations  (8)  and  (9)  are  the  equations  of  transfer .  The  boundary 
conditions  for  these  equations  depend  on  the  particular  problem  at  hand. 

In  the  case  of  an  absorbing  and  emitting  medium,  equations  (8)  and  (9)  are 
simplified  considerably. 

In  such  a  case 

x 

6A  “  °A  and 


9 


and  (8)  and  (9)  become 


w 


(12a) 


(12b) 


The  monochromatic  radiation  in  the  X-direction  is  given  by 

V(tx>  '  /  l,d" 

4* 


(13) 


for  a  one  dimensional  radiative  transfer  case,  the  boundary  conditions  for 
equations  (8),  (9)  and  equations  (12)  can  be  written  as: 


«  ta  -  0 


at 


(14) 


The  expression  for  0  for  absorbing,  emitting  and  scattering  media  becomes 

R 

qm(’x)  "  2' J  / I‘('01.-i.)e'(,0x',x)/l,u<lu 


* 2  ♦srGx(t>]  ‘2<vt)dt 
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- 2  /  [irebx(t)  *4S7-Gx(t)]^(t-’.)dt 
1 

X 


(15) 
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and 

/in  ^  1  ' 

I'tO.wJe-V"  d„  ♦  2,f 

0 

(18) 

r'0x 

*2'J  *b»(t)  «i<l'A-t|)dt  -4  *bl(t) 

Equations  (8),  (9).  and  (14)  are  quite  difficult  to  solve  numerically  and 
result  in  a  very  complex  computer  program. 

An  alternative  analysis  has  been  developed  using  a  different  approach 
which  results  in  considerable  simplification  in  the  radiative  transfer 
analysis  for  an  absorbing,  emitting,  and  scattering  medium.  This  develop¬ 
ment  is  described  in  the  next  section. 


1 1 


FIGURE  3  COORDINATE  SCHEME  FOR  RADIATIVE  HEAT  TRANSFER  BETWEEN 

TWO  INFINITE  PARALLEL  PLATES 
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B.  Ef fectlve  Absorptivity  Anol ysls 

It  con  be  seen  from  the  previous  section  thot  by  ossumlng  the 
scotterlng  coefficient  y  to  be  zero,  and  therefore  .  the  equations 

8.  9,  and  15  can  be  simplified  considerably  to  become  equations  12a,  12b 
ond  17  respectively.  The  simplified  equations  12  and  17  neglect 
scotterlng.  If,  however,  an  expression  for  the  obsorption  coefficient 
can  be  found  thot  Includes  the  effect  of  scattering,  then  equations  12  ond 
17  will  describe  radiative  heat  transfer  in  an  absorbing,  emitting  and 
scattering  medium.  This  is  precisely  what  Is  done  in  the  Effective  Absorp 
tlvlty  Analysis  described  below: 

Consider  a  case  of  radiative  heat  transfer  between  two  infinite 

parallel  plates  with  an  intervening  medium  of  on  absorbing,  emitting  and 

scattering  gas  (Figure  3).  Consider  a  thin  layer  of  the  gas  in  the  medium 

Assume  the  radiant  heat  flux  vector  to  be  In  the  form  of  a  difference  of 

two  fluxes  in  the  X-dlrection  flowing  in  opposite  directions.  The 

♦ 

differential  equations  for  the  Intensity  of  rodiotion  ond  1^,  as 
derived  earlier  can  be  written  as: 

*  w~  ~  ■  *  sr/  1x<M'’x)s<‘*'-u)d"'  * 


dl~(p»x) 


-6xI~(u,x)  ♦  Ix(jj',x)S(u',u)dw' 


axebx 


Integration  of  (B)  over  a  hemisphere  gives: 


dw  »  -a^  J  Ix(u,x)dw  -  y yj  I A ( ij . x)dw 
2%  2% 


_iL .Cdutji '  ,x)S(u'  ,u)du>'  ♦  2t>xe^ 


Here  I  represents  integration  over  hemispherical  solid  angle  (Figure  3) 
*'2*4 

sucn  that  g.  measured  from  the  positive  normal,  varies  from  1  to  0  and  / 
represents  integration  over  a  hemispherical  solid  angle  such  that p  measured 
from  the  positive  normal  varies  from  0  to  -1. 

The  second  term  on  the  right  hand  side  of  (19)  con  be  written  as 


"T»  /  II(u>x)d"  ■  f 


2*. 


'4»  A<->  / 

V  4* 


f S(p  ,p ')dw  1  dw' 


[  ^  f  S(p.p ')dw=l  J 


s  ~  f  dw  y*  IA(u',x)S(p,p')dw' 
4s  2«+ 


(20) 


Substituting  (10)  in  (19).  we  get  the  right  hand  side  equal  to 
-aX  1*  (w.x)dw  -  j " dw  /”lx(u',x)S(p,y ')dw 


2%  2% 


+  /d"/  IA(w'.x)S(p,p')dw'l+  [/  dm  1^  (v  '  ,x)S(p  ,p  ')dui' 

2s_  2w+ 

*  J  dw  /  I  (p',x)S(p,p')dw'”|  +  2a 

•V  y2s  X  J  A 


2*+  2*+ 


Noting  that  integration  of  1^  on  2X  will  be  the  same  as  integration  of 
♦  ♦  + 
on  2  ,  equation  (19)  becomes 

5x  /  f  >»<».*)«<»  -  JV /  d“/ (J1) 

2%  2*+  2*_  "Zir  + 


£/*■/ 


2*+ 


I .  ( u ' ,x )S ( ii , ii '  )dw '  ♦  2a  e. 

*  a  dx 


14 


Equation  (21)  can  be  written  in  compact  form  as: 


VE».*r 

Similarly,  from  (8b)  one  can  get: 


(22a) 


ir  *  K  *  V»  V». '  V»V»*  * 2 


(22b) 


q 


Equations  (22o)  and  (22b)  are  the  differential  equations  for  the  fluxes 
E  and  E  In  these  equations,  the  terms  are  defined  as  follows: 


X  +  X. 


E».  *  /  ■» 

t  «  f  lx  Md» 


Spectral  volumetric  density  of  spontaneous  radiation. 


I 

/»■ 

.1 


I 

ft1 

»' 


,<  j 
<1 
•Y 

rf  f 


4<XX°bX 


J 


/  I  (v,x)gdw 

'2*+  X 

^  I^(p,x)dw 

2w_ _ 

I^(u,x)gdw 

2* 

/  d ul  I  (u',x)S(g,g')d<, 

J2r  J2*a  X 


4n  /  I  (w'.x)dw' 
'2%  A 


/  <k>  /  l ,  (  M ' ,x)S(p ,p ')dw' 

'2%  y2i>_  A 

4«  /  I  (p ',x)dw' 

/2«  A 


If  I  Is  isotropic  within  the  limits  of  the  hemispherical  angle  m  and  m 
X  A  +  X_ 

turn  out  to  be  equal  [8j. 

m  *  m  =  m  =  2 
A+  A_  A 

while  the  coefficients  6  and  6  for  the  axisymmetrical  characteristic 

+ 

curves  of  arbitrary  shape  are  also  found  to  be  equal  to  one  another  and  are 


determined  from  the  following: 


6.  -  6.  *  6  *  - 

x-  X  8* 2 


—  /  dw  / 
fl.2  J  ♦  l 


S(w,y  ')dw' 


2w  2s 


From  (22a)  and  (22b) 


d(V  V 


-2(y2V»l£i  *2<«x+2Vx>Ex  *  R» 


*  -*V2VxK  *  Rx 


dcT 

- -  .  —  •  -2axE, 


-*vl 


ET  *  "2(aX  *  2ViK  *  RX 


2b  ■  -2axcx 


from  (31)  and  (32)  one  obtains 


«xS  *  ‘2oxRx 


2“xEx. 


iV* 


:,nv 


V  * 


where 


'  °j<0j  *  2Vi> 


The  solution  of  this  equation  is  found  to  be 


2K  x  -2K  x  a  R. 
c,  =  Ae  A  +  Be  A  ♦  - 


and 


(34) 


(35) 


+ 


] 


Therefore 


<,"«  ea 


are  obtained  as 


(36) 


(37) 


(38) 


Boundary  conditions  for  the  layer  are 


and 


Substitution  of  these  conditions  in  (37)  and  (38)  gives  A  and  B  as 


of  the  layer  is 


Reflectivity  "R  " 
f 


Absorptivity  “A 

SC 


D  -  R 


f 


(45) 


Equation  (45)  is  the  expression  for  absorptivity  which  includes  the  effect 
of  scattering.  This  is  given  the  name  “effective  absorptivity  in  the 
presence  of  scattering . " 

Recalling  that 

Absorptivity  A(x)  *  (46) 


an  effective  absorption  coeff iclent  a.  can  be  obtained  from  the  effective 


sc 


absorptivity  of  a  layer  of  thickness  (x  -x  )  from  the  relation 

2  1 


\c  ■  '  -  c 


I’ 


(47) 


a.  thus  calculated  is  the  coefficient  which  takes  into  occount  the  effect 
*sc 


of  scattering . 


Absorption  coefficients  calculated  by  this  method  ore  used 


In  the  calculation  of  heat  radiation  from  equations  (17)  and  (18): 


and 


ydp 


1  /„  1 
=  J"lx(o,\i)e  *  wdy  -  2 nf _y)e 

0  *0  (17) 

Tx  T°A 

*  2f  -  2 /  ebx(t>c2(t-^)dt 


dQ 


..  2,  /"  *  2 


I 


/OX  /  \ 

•bx(t)M(|,x*t|)  -  4ebx(’x> 


(18) 


in  these  equations  is 


“xscdx 


(US) 


It  should  be  noted  that  the  above  equations  (17  and  (18)  are  for  an 
absorbing  and  emitting  medium.  By  substituting  the  coefficient  a 

sc 

calculated  from  (45)  and  (47)  for  the  regular  absorption  coefficient  , 

these  equations  will  describe  the  radiation  transfer  for  an  absorbing, 
emitting,  and  scattering  medium. 

Assume  that  the  wall  surface  is  black  with  respect  to  radiation,  and 
that  any  radiation  incident  on  the  edge  of  the  boundary  layer  has  a 
wavelength  distribution  identical  to  a  perfect  radiator,  and  furthermore, 
that  radiation  emitted  from  the  surface  or  at  the  edge  of  the  boundary 
layer  is  emitted  diffusely.  Then 
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4 


Scotterlnq  Function 


The  scottering  function  (  6 ,  0  ;0  ,♦  ),  sometimes  colled  phose 
function,  characterizes  the  directional  distribution  of  the  scattered 
energy,  such  that 

sx(e',4>';  6,  ♦  )  47 

represents  the  probability  that  radiation  in  the  frequency  interval  d^ 
centered  ot A  and  from  a  direction  (0'$')  will  be  deflected  in  a  direction 
(0  ,  9  )  within  a  solid  angle  "q^".  Figure  4  shows  a  coordinate  scheme  for 
the  scattering  function.  The  integral  of  this  quantity  over  all  the  solid 
angles  should,  therefore,  be  equal  to  1. 

i.e.  Sx(e',*';e,*)  d«  =  1 

4* 

Scattering  by  a  single  particle  is  called  "single  scottering."  If  a 
relatively  few  particles  exist  in  the  volume,  the  scattering  can  still  be 
characterized  by  "single  scattering"  but  will  have  an  amplitude  proportional 
to  the  number  of  particles  in  the  volume.  As  the  number  of  particles  is 
Increased,  each  particle  is  subjected  to  increased  radiation  scattered  by 
other  particles.  This  type  of  scattering  is  characterized  by  "multiple 
scottering . " 

Analytical  determination  of  the  scattering  function  is  made  by  solving 
Maxwell's  Equations  for  the  propagation  of  electromagnetic  waves  with  the 
appropriate  boundary  conditions.  The  problem  was  first  solved  by  Gustav 
Mie.  The  development  and  results  were  reviewed  in  detail  by  Van  De  Hulst 
(0)  ond  will  not  be  given  here.  According  to  Mie's  theory,  for  homogeneous 
spheres,  the  scottering  function  may  be  expressed  as  a  function  of  the  size 
parameter  "2"  ond  the  refractive  index  "m-ix”  of  the  sphere  with  respect  to 


22 


FIGURE  4.  COORDINATE  SCHEME  FOR  SCATTERING  FUNCTION 


the  surrounding  medium. 


Z  = 


(55) 


where 

D  ■  Diometer  of  the  particle 

P 

A  -  Wavelength  of  radiation 

For  values  of  2«1,  a  simpler  formulation  developed  by  Rayleigh  holds 
and  for  large  values  of  2,  principles  of  geometrical  optics  ore  often 
utilized.  The  computation  of  the  scattering  function  is  thus  frequently 
classified  as  below: 

2  «  1  Rayleigh  scattering 

2  a  1  Mie  scattering 

Z  »  1  Geometrical  scattering 

When  the  scattered  radiation  is  spread  equally  in  all  directions,  it  is 
known  as  isotropic  scattering. 

Sx(e',  0,  f)  *  I  (56) 

The  scattering  function  for  RAyleigh  scattering  can  be  written  in  the 

form 

sx(e',*';e,*)  *  l+Cos2e  )  (57) 


Here  0  represents  the  polar  angle  between  the  incident  ray  and  the  scattered 
ray. 

In  this  study,  axially  symmetric  scattering  functions  have  been  used. 

An  axially  symmetric  scattering  function  is  defined  as  describing  the 
condition  in  which  the  intensity  of  radiation  scattered  from  a  pencil  of 
rays  incident  on  an  elemental  volume  is  axiolly  symmetric  with  respect  to 
the  direction  of  the  incident  pencil.  A  further  requirement  is  that  the 
scattered  intensity  be  independent  of  the  direction  of  the  incident  pencil 


with  respect  to  the  axis  of  the  system.  Thus  for  the  case  of  o  particle 
with  an  axially  symmetric  scattering  function,  the  function  may  be 
tabulated  os  a  function  of  o  single  angle  (usually  the  angle  between  the 
scattered  ray  and  the  forward  direction  of  the  incident  pencil  of  rays). 

The  assumption  of  axially  symmetric  scattering  function  has  been  shown  to 
be  a  reasonable  assumption  for  random  orientation  of  particles  by  Love  and 
Beattie  [9].  They  presented  measurements  for  clouds  of  aluminum  particles, 
iron,  glass,  and  silica  particles.  A  typical  plot  of  the  experimentally 
obtained  scattering  function  is  given  in  Figure  5. 

In  this  study,  the  work  of  Chu,  Churchill  and  Clark  [10]  was  used  to 
compute  the  scattering  functions.  They  prepared  their  values  for  real 
refractive  index  and  axially  symmetric  scattering.  The  extinction  and 
scattering  coefficients  were  taken  from  the  work  of  Chromey  [ll]  and  are 
presented  in  Figures  6  and  7. 

Chu.  Churchill  and  Clark  [10]  have  published  a  table  of  coefficients 

to  be  used  in  a  Legendre  Polynomial  representation  of  the  scattering 

function  for  spheres  having  real  refractive  Indices  ranging  from  m  -  0.9  to 

m  -  2.0  and  7R  »  «>  and  for  values  of  Z  from  1  to  30.  In  their 

representation,  the  scattering  function  is  expressed  as 

n=« 

s(0)  =  1  +  X  an  Pn(Coso)  (58) 

where  S( 0)  is  the  axially  symmetric  scattering  function,  o  are  tabulated 

n 

coefficients  dependent  on  the  particle  size  parameter  "Z"  and  the 

refractive  index  "m."  p  (Cos0)  are  the  nth  order  Legendre  Polynomials  with 

n 

the  argument  Cos9  and  0  is  a  polar  angle. 
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SCATTERING  FUNCTIOI 


POLAR  ANGLE  6 


FIGURE  5.  A  TYPICAL  PLOT  OF  ANGULAR  DISTRIBUTION  OF 


SCATTERING  FUNCTION  [7] 


FIGURE  6.  THE  PARAMETER  YX/0X  AS  A  FUNCTION  OF  PARTICLE  SIZE  PARAMETER  Z 


Rodlotlve  Properties  of  Surfoces  4  Portlcles 


Radiative  properties  that  hove  importance  in  the  high  speed  vehicle 
development  include: 

1.  Emittonce  and  reflectance  of  vehicle  skin  surface 

2.  Absorption  coefficient  and  scattering  coefficient  of 
particles  that  can  be  injected  into  the  boundary  layer 
as  coolants. 

3.  Absorption  and  scattering  coefficients  of  ablation  products 
of  vehicle  sheathing 

4.  Scattering  functions  of  coolant  particles  ond  products  of 
ablation . 

A  computer  literature  search  was  conducted  to  compile  the  published 
data  of  radiative  properties  of  materials  important  to  this  study.  The 
available  data  is  being  presented  in  this  report.  There  is  very  little 
data  available  on  the  radiative  properties  of  materials  of  interest.  Most 
of  the  data  available  is  for  the  emisslvlty  and  reflectivity  of  surfaces  of 
metals  [1-93.  The  data  for  absorption  and  scattering  properties  of 
particles  is  almost  nonexistent.  Goubareff,  et.  al.  [12]  compiled  the 
thermal  radiation  properties  data  available  in  the  fifties  ond  before.  The 
data  reported  was  mostly  for  metals,  alloys  and  building  materials. 
Touluklan  and  DeWitt  [13]  also  compiled  the  radiative  properties  of  metals 
and  alloys.  Unfortunately,  most  of  the  data  reported  in  References  1  and  2 
is  old  and  there  is  considerable  disparity  between  data  token  by  different 
investigators  for  the  same  materials.  Blau  and  Francis  [14]  presented 
measured  data  of  the  spectral  emittonce  of  stainless  steel  and  platinum  and 
inconel  surfaces  coated  with  silicon  monoxide  for  a  wavelength  range  of  2 
to  14  microns.  They  also  showed  the  effects  of  surface  oxidation  which  con 
change  the  spectral  properties  by  as  much  as  30M.  The  data  is  presented 
in  Appendix  A.  Neuer  and  Warner  [15]  also  showed  the  effect  of  surface 
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oxidation  for  steel.  Keegan,  et.  ol.  [16]  presented  measured  data  on  the 
diffuse  spectral  reflectance  of  beryllium,  steel,  aluminum  and  platinum 
over  a  wavelength  range  of  .2  to  2.1  microns.  This  data  is  also  presented 
in  Appendix  A.  PeletskiJ.  et.  al .  [17]  presented  data  on  the  total  hemispheri¬ 
cal  emlttance  of  titanium.  The  more  recent  studies  on  radiative  properties 
of  metallic  materials  are  presented  in  References  18  and  19.  Ramanathon  [18] 
compared  some  recent  measurements  with  theoretical  predictions.  The  data 
for  copper,  silver,  and  aluminum  is  presented  in  the  Appendix  A.  Makiio, 
et  al.  [19]  evaluated  theoretically  the  spectral  and  total  emisslvitles  of 
cobalt,  chromium  and  nickel.  Their  calculated  values  for  cobalt  compared 
well  with  the  measured  values.  Therefore,  their  data  for  cobalt,  nickel 
and  chromium  is  presented  In  the  Appendix  A. 

Radiative  properties  of  surfaces  of  ceramics  and  other  heat  resisting 
materials  are  presented  in  References  20  to  27.  DeWitt  and  associates  [20. 

21]  presented  measured  values  of  spectral  emissivlty  of  silicon  carbide, 
silicon  nitride  and  tantallum.  These  values  are  given  in  the  Appendix  A. 

Wilson  [22]  reported  data  of  spectral  and  total  reflectivity  and  emissivlty 
of  ablation  chars  and  carbon  and  zirconla  for  temperatures  of  2100 'K  to 
3700 *K.  This  data  is  also  presented  in  Appendix  A.  Douglas  [23]  has 
presented  the  radiative  properties  of  a  Mylar-aluminum  laminate.  Khrustaler 
and  Rakor  [24,  25]  presented  measured  data  for  total  hemispherical  and 
total  directional  emissivlty  for  Noblum,  Molybdenum  and  tontalum.  The 
most  recent  dota  of  radiative  properties  of  ceramic  materials  was  presented 
by  Maklno,  et.  al .  [26,  27,  28].  They  presented  spectral  and  total 
properties  of  surfaces  of  white  ceramics  (Al  0  ,  Z  0  ),  black  ceramics  (Si 
N  ,  SIC)  and  metallic  ceramics  (TIC  and  TIN).  Their  dota  as  well  as  data 
from  references  23  and  24  is  presented  in  the  Appendix  A. 
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Shofey  and  Kunltomo  [29,  30]  and  Kunltomo,  et.  al .  [31]  have  presented 
theoretical  and  experimental  studies  of  the  absorption  and  scattering 
properties  of  particles  of  pigments  and  reloted  these  properties  to  the 
reflectance  of  painted  layers  containing  those  pigments.  While  the  data 
presented  is  not  of  much  importance  to  our  study,  some  of  the  ideas 
presented  in  thel’*  papers  are. 

The  data  that  is  lacking  the  most  in  the  literature  is  for  the 
obsorptlon  and  scattering  properties  of  particles  and  clouds  of  particles. 
References  32  and  33  presented  the  radiative  properties  of  clouds  of  soot 
generated  by  fires.  While  this  data  is  not  directly  useful  in  this  study, 
it  may  become  useful  if  a  method  of  thermal  protection  is  developed  based 
on  this.  Therefore,  selected  data  from  these  references  is  given  in  the 
Appendix  A.  Reference  34  to  36  presented  some  useful  data  on  the  radiative 
properties  of  particles.  Nagy  and  Lenoir  [34]  presented  absorption  and 
scattering  coefficients  of  aluminum  oxide  particles  of  sizes  6  and  12 
microns  In  the  wavelength  range  of  2  to  12  microns.  Sanders  and  Lenoir 
[35]  presented  measurements  of  back  scattering  ond  extinction  efficiency  of 
graphite  and  aluminum  oxide  particles.  Love  ond  Beattie  [36]  also 
presented  measured  data  for  clouds  of  aluminum,  Iron,  glass  and  silica 
particles.  Some  selected  data  from  these  references  In  presented  in  the 
Appendix  A. 


RESULTS 


Based  on  the  theoretical  development  described  earlier  in  this  report, 
o  computer  program  was  developed  that  solves  the  turbulent  boundary  layer 
flow  equations  for  mass  injection  of  the  obsorping,  emitting  and  scattering 
particles  in  the  boundary  layer.  This  was  done  by  modifying  an  earlier 
code  developed  by  the  author.  This  earlier  code  [l]  involved  the  inclusion 
of  radiative  heat  transfer  analysis  in  the  Patankar  Spalding  code  for 
turbulent  boundary  layer  flow.  The  modification  involved  the  Inclusion  of 
concentration  of  species  equation  in  the  boundary  layer  flow  to  account  for 
the  mass  injection  and  ablation.  This  modification  was  taken  from  another 
publication  by  the  author  (Reference  3).  Since  both  of  the  above 
mentioned  studies  (l.  3)  had  validated  their  results  separately,  and  the 
present  study  involved  merging  the  two  no  attempt  was  made  in  comparing  the 
present  results  with  any  other  publication. 

The  developed  computer  program  was  used  to  obtain  results  for 
radiative,  convective  and  total  heat  transfer  to  the  surface  of  a  flat  wall 
for  various  combinations  of  wall  temperature,  free  stream  temperature,  Mach 
Number  and  scattering  coefficients  of  the  injected  particles.  The  results 
were  obtained  to  analyze  the  effect  of  injected  particles  for  protecting 
the  wall  surface  from  incident  thermal  radiation  from  an  external  source. 

Figure  8  shows  the  total  heat  flux  in  Btu/ftl-sec.  to  the  surface  as  a 
function  of  the  distance  from  the  leading  edge  of  the  plate  for  T^/T  of  1, 
2,  4  and  6.  The  total  heat  transfer  is  the  sum  of  the  convective  and 
radiative  heat  transfer.  This  figure  shows  the  totol  heat  transfer  when 
there  is  no  external  incident  radiation.  In  such  a  case,  the  radiative 
heat  transfer  is  a  very  small  part  of  the  total  heat  transfer.  This  con  be 
seen  from  Figure  9.  As  seen  from  this  figure  the  radiative  heat  transfer 
is  less  than  .  "\%  of  the  total  heat  transfer  for  T^/T^  equal  to  1.  The 


Fig:  8  TOTAL  HEAT  TRANSFER  AS  A  FUNCTION  OF 
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DISTANCE  FROM  LEADING  EDGE  (ft. 


Fig:  9  RADIATIVE  AND  TOTAL  HEAT  TRANSFER  AS  A  FUNCTION  OF  FREE  STREAM  TEMPERATURE 
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relative  fraction  of  the  radiative  heat  transfer  Increases  os  T  Increases 
for  the  some  .  At  T^/T^  of  6  the  percentage  of  the  radiative  heat 
transfer  Increases  to  about  5if.  Figure  10  shows  the  radiative  heat  flux  to 
the  surface  os  a  function  of  the  distance  from  the  leading  edge  for  various 
free  stream  temperatures .  It  is  seen  from  these  two  figures  that  radiative 
heat  flux  is  relatively  unimportant  when  the  free  stream  temperatures  ore 
not  too  high.  In  such  a  situation  any  strategy  for  thermal  protection 
should  try  to  reduce  the  convective  heat  flux.  However,  protection  from 
thermal  radiation  becomes  important  when  there  is  a  strong  external 
incident  radiation  source  present,  such  as  a  nuclear  blost  or  a  laser  beam. 

In  this  study,  we  have  assumed  the  external  radiation  source  to  cause  a 
uniform  incident  thermal  radiation  flux  of  1000  Btu/ft*-sec.  at  the 
external  edge  of  the  boundary  layer  (Figure  11).  Figure  12  shows  the 
effect  of  Mach  Number  on  the  radiative,  convective  and  total  heat  flux.  As 
seen  from  this  figure,  the  radiative  heat  flux  is  affected  very  little  by 
the  change  in  Mach  Number.  However,  the  convective  heat  flux  increases 
with  Mach  Number,  thereby  increasing  the  total  heat  flux.  Figure  13  shows 
a  typical  temperature  profile  in  the  boundary  layer  in  hypersonic  flow. 

The  injection  of  particles  that  absorb  and  scatter  thermal  radiation 
reduces  the  externally  Incident  radiative  flux.  This  is  seen  from  Figures 
14  and  15.  In  these  figures  the  effect  of  scattering  coefficient  of  the 
injected  particles  on  the  radiative  and  total  heat  flux  is  seen.  These 
figures  show  very  clearly  the  potential  benefits  of  using  this  technique 
for  the  thermal  protection  of  the  vehicle  surfoce  by  injection  of  strongly 
scattering  particles.  Figure  16  shows  the  effect  of  scattering  on  the 
radiative  heat  transfer  for  different  Mach  Numbers.  It  can  be  seen  from  this 
figure  that  by  injection  of  purely  scattering  porticles  it  is  possible  to 


Fig:  10  RADIATIVE  HEAT  TRANSFER  AS  A  FUNCTION  OF  FREE  STREAM  TEMPERATURE 
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Fig:  12  RADIATIVE  AND  TOTAL  HEAT  TRANSFER  AS  A  FUNCTION  OF  MACH  # 
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Fig:  13  A  TYPICAL  TEMPERATURE  PROFILE  IN  THE  BOUNDARY  LAYER 
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SCATTERING  COEFFICIENT 


Fig:  15  EFFECT  OF  SCATTERING  COEFFICIENT  ON  RADIATIVE  AND  TOTAL  HEAT  TRANSFER 


SCATTERING  COEFF 


MACH 


reduce  the  Incident  radiative  flux  by  75^  for  flow  at  Mach  Number  of  6. 
Similar  reductions  are  seen  at  higher  Mach  Numbers  also.  Figure  17  shows  o 

temperature  profile  in  the  boundary  layer  with  the  injection  of  absorbing 
and  scattering  particles.  If  the  particles  only  absorb  and  do  not  scatter, 
the  temperature  in  the  boundary  layer  will  be  higher  than  the  case  where 
particles  scatter  radiation.  The  reason  is  that  for  the  case  of  scattering 
particles  part  of  the  incident  radiation  is  reflected  out  of  the  boundary 
layer  instead  of  being  absorbed  resulting  in  lower  temperatures  in  the 
boundary  layer.  The  results  show  that  for  the  case  where  there  Is  no 
external  incident  thermal  radiation  the  radiative  flux  is  a  very  small 
fraction  of  the  total  heat  flux.  In  this  case,  the  injection  of  mass  does 
not  have  any  effect  on  the  radiative  heat  transfer.  However,  even  in  such 
a  cose,  there  is  a  benefit  that  can  be  seen  from  Figure  18.  This  figure 
shows  that  mass  injection  in  this  case  reduces  the  convective  and  thus  the 
total  heat  flux.  This  happens  because  the  injection  of  mass  reduces  the 
temperature  gradient  of  the  wall,  thereby  reducing  the  convective  heat 
flux . 
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APPENDIX  A 

RADIATIVE  PROPERTIES  OF  METALS,  ALLOYS,  CERAMICS,  ABLATION  CHARS, 
SOOT  CLOUDS  AND  PARTICLES,  COMPILED  FROM  PUBLISHED  LITERATURE 

NOTE:  Figure  captions  are  taken  from  original  papers.  Reference 
numbers  and  equation  numbers  in  figures  refer  to  those 
in  the  original  papers. 
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SPECTRAL  EMITTANCE  OF  STAINLESS  STEEL,  PLATINUM,  INCONEL 

(Data  from  Ref.  14) 


Finiinr.  2— Spectral  tmiUanee  •/  eilican  mcnozitU  conted 
platinum  and  Inconel. 


REFLECTANCE  OF  BERYLLIUM,  STEEL,  ALUMINUM  AND  PLATINUM 

(Data  from  Ref.  15) 


Fioomc  1. — Diffuie  tpedral  reflectance*  ( 0.16  to  t.l  «■) 
of  machine-lapped  and  eandblaeted  tample*  of  alu¬ 
minum  and  of  platinum.  Back  euro*  repmenl* 
determination*  at  four  orientation*.  Tkert  too*  no 
appreciable  effect  of  orientation. 


of  eurfaoc  of  beryllium  eample  1,  bond-lapped  to  about 
4  microinche*.  The  tune  represent*  determination*  at 
eight  orientation*  of  the  t ample. 


Fiounc  5. — Diffute  tpectral  reflectance*  (O.tt  to  t.l  *) 
of  iteel  gage  block*  with  commercial  flniihe*.  The 
variation  in  reflectance  it  due  to  difference*  in  the 
•tevle  need  in  the  manufacture  of  the  block*  and  to 
a*  id*  turf  ace  film*. 

I.  Webber  !i  Min  St-t  l.  Fan  Keuren  0.150  in. 
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SPECTRAL  EMISSIVITY  OF  SILICON  CARBIDE  AND  SILICON  NITRIDE 

(Data  from  Ref.  20,  21) 


rig.  3  Homl  spectral  oduivlty  of  •illeon 
carbide  at  1900K 


rig.  S  Normal  spectral  aolsslvlty  of  silicon 
nitride  at  1700K 
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SPECTRAL  AND  TOTAL  EMITTANCE  OF  PHENOLIC  NYLON  CHAR 
(Data  from  Ref.  22) 
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SPECTRAL  AND  TOTAL  EMISSIVITY  OF  MOLYBDENUM  AND  NOBIUM 
(Data  from  Ref.  24) 


e 

a 

KUO  MO  >too  iuo\ 


A 

A 

of 

aJ 

-&>***} 

-v — 

Fig.  1.  The  emissivlties  of  molybdenum  as 
a  function  of  temperature,  a)  Total  hemi¬ 
spherical  emisalvlty;  b)  total  normal  emts- 
sivity;  c)  spectral  emisalvlty  for  A  *  0.65  n; 
1)  experiments  In  vacuum;  2)  experiments 
In  argon;  3)  preheating. 
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Fig.  2.  The  spectral  cmissivity  of  molybdenum  as  a  func¬ 
tion  of  the  wavelength.  1)  Measured  by  infrared  spectrom¬ 
eter  at  t  -  1681*C;  measured  by  infrared  spectrometer  at 
t  =  1825*C;  3)  measured  by  the  photographic  method;  4) 
measured  by  the  optical  pyrometer. 


Fig.  3.  Integral  emlssivity  of  a  niobium 

specimen  (tubes  with  d  .  »  6. 29.  d.  * 
oui  in 

5. 29)  for  studies  in  vacuum  and  argon.  1) 
Studies  In  vacuum;  2)  studies  in  argon  (with¬ 
out  correction  for  natural  convection) ;  3) 
studies  in  argon  (with  correction  for  natural 
convection). 
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SPECTRAL  AND  TOTAL  EMISSIVITY  OF  TANTALUM 


(Data  from  Ref. 


21  ) 


(Data  from  Ref.  24) 


Fig.  7.  The  emlssivities  of  tantalum  as  a 
function  of  temperature,  a)  Total  hemispher¬ 
ical  emt88(vtty;  b)  total  normal  emtssivity;  c) 
spectral  emissivlty  for  A  =  0. 65  ;  1)  experi¬ 
ments  under  vacuum;  2)  experiments  in  argon; 

3)  Initial  heating. 


ABSORPTION  COEFFICIENT  AND  EMISSIVITY  OF  SOOT  CLOUD 
(Data  from  Ref.  32) 


WAVELENGTH  A  ym 


FIG.  3  Spectral  absorption  coefficient  of  a  cloud  of  small  soot 
particles. 
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ABSORPTION  AND  SCATTERING  PROPERTIES  OF 
ALUMINUM  OXIDE  PARTICLES 
(Data  from  Ref.  34) 
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APPENDIX  B 

COMPUTER  PROGRAM  LISTING 


c 


MAIN  PART  OF  THE  PROGRAM 

C0MM0N/BLK2/N »  NP1 »  NP2 f  NP3  f  NEQ  f NPH »  KEX f  KIN  f KASE  f  KRAD »  PE I »  AM I r  AME  » 

1  DPDXfXUfXDfXPfXL  fDXf INTGfCSALFAfHGSfCINF< 10) fPREF(IO) » 

2  PR< 10) fP(10) »  U ( 050 ) , F( 10 1 050) fR'050) fRH0(050) f0M<050) f Y(050) » 

3  DEN(050) »AMU(050) 

C0MM0N/BLK3/KRPfMACHfKRf INDI ( 1 0 ) » I NDE < 10 ) ?  SUMR  t AK »  ALMG f BETA  f  TAU I 

1  fTAUEfVINFfDELTAfTWfTINFfGSfYLfUMAXfUMINfFRf YIPfYEMfTRfHWf 

2  HF< 10) t GAMA  < 10)fAJI<10)fAJE(10)fSC(050)fAU(050)f  BU(050) fCU(050) r 

3  A ( 10f050) f  B  < 10r050) »C< 10»050) »TEMPE(050) f  TEMP ( 050 ) f  PO ( 050 ) » 

4  AMACH(050) fCP<050) 

C0MM0N/BLK4/NHfNN0fNNPfN0fNNfN02fNN2fNEfIND( 10) »  STO  »  AKS  »  RT fFTf 
1  AMTfDRAGCfSTNOfTAUIUfOUfUGUfUGD 

C0MM0N/BLK5/KRI f  PINF ,DIFF (50) fNOPTfNSKAT 
C0MM0N/GC/GC0N<050> fAMOL (050) 

CHARACTER*16  NAME 

TYPE#  t  •  TYPE  THE  NAME  OF  THE  OUTPUT  FILE' 

ACCEPT201 fNAME 
201  FORMAT  <  A 16 ) 

0PEN<UNIT  =  1fFILE=NAMEfSTATUS=/NEW'  ) 

N  =  30 

nopt  -  0 

10  CONTINUE 

X  =  0.0 
I NTG  =  0 
CALL  CONST 
GO  TO  20 
15  N  a  N+l 

INTG  =  0 
XU  =  0.03 

KR  =  KRI 

20  CALL  BEGIN 

CALL  EXTRA 
AM  I  =  0 . 

AME  ~0 . 

GO  TO  30 

25  CALL  READY 

30  CONTINUE 

I  NTG  =  I  NTG-f  1 
IF  < INTG . GT . 50 ) KR=KRI 
CALL  LENGT 
CALL  ENTRN 

C  CHOICE  OF  FORWARD  STEP 
FRA=0 . 05 

DX=FRA*PEI/(R( 1 )*AMI-R(NP3)*AME ) 

I F ( D  X . G  T . . 5* Y ( NP3 ) > DX= . 5*Y < NP3 > 

XD=XU+DX 

C  ASSUMES  FREE  STREAM  VELOCITY  =  AX  +  BXX  ♦  C 

C  CALCULATES  CHANGE  IN  FREE  STREAM  VELOCITY 

A 1  =  0.0 
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Bl  =  0.0 
Cl  =  VINE 
U6U  =  UGP 

UGD  =  A1  *  XD  +  B 1  ♦  X  D  *  XD  +  Cl 

DPDX=  < (JGU  +  UGD  )  *  (  UGU-UGD  >  *  .  5*RH0(NP3  )  /  (XD-XU) 

IF  < KASE . EQ .  2 )  GO  TO  35 

IF (KIN. EG. 1 ) CALL  MASS < XU . XD » AM I ) 

IFIKEX.EQ. 11CALL  MASS < XU . XD » AME ) 

CALI.  WALL 
35  CALL  COEFF 
CALL  OUTPUT 

C  SETTING  UP  VELOCITIES  AT  A  FREE  BOUNDARY 
IF  (KEX.EQ.2)  U < NF‘3 )  =UGD 

IF(KIN.EQ.2)  U( 1 )=SQRT(U( 1 ) -2 . * ( XD-XU ) «DPDX/RHO < 1 ) > 

CALL  SOLVE<AUf BUf CU.U»NP3) 

C  SETTING  UP  VELOCITIES  AT  A  SYMMETRY  LINE 
IF<KIN.NE.3>  GO  TO  40 
U< 1 )*U<2> 

IF  <  KRAP . EG . 0 )  U< 1 >  =  . 75*U < 2 > + . 25*U < 3 > 

40  IF  <  KEX . EQ . 3  )  U < NP3 > = . 75*U < NP2 > + . 25*U ( NPl > 

45  CONTINUE 

IF(NFG.EQ.l)  GO  TO  90 
DO  75  J=  1 f  NPH 
DO  50  I=2»NP2 
AU< I >*A< J. I  ) 

BU( I)=B<J?I) 

50  CU<I)=C<J.I> 

DO  55  1  =  1 f  NP3 
55  SC<I)=F(J»I) 

CALL  SOLVE (AUfPUrCUrSC?  NP3  > 

DO  60  1=1 »NP3 
60  F(J,I)«SC(I) 

IF  <  KASE . EQ . 2 )  GO  TO  65 
C  SETTING  UP  WALL  VALUES  OF  F 

IF(KIN.EQ.1.AND.INDI(J).EQ.2)F(J«I)=(( 1 . +BETA+GAMA ( J ) )*F( J,2) 
1  < 1 . 4 BET A -GAMA ( J> )*F ( J»  3  ) )*. 5/GAMA (J) 

IF (KEX. EG. 1 .AND. INDE< J) . EQ.2)F( J,NP3>*< ( 1 . +BETA+ GAMA  (  J ) )* 

1  F( J»NP2)-(1 ,+BETA-GAMA( J) )*F( J»NP1 >>*. 5/GAMA < J ) 

C  SETTING  UP  SYMMETRY-LINE  VALUES  OF  F 
65  IF<KIN.NE.3>  GO  TO  70 

IF(KRAD.EQ.O)F< J»1 >«.75*F( J.2>F.25*F( J. 3) 

70  IF  <  KEX . EQ . 3  >  F ( J . NP3 ) * . 75*F  <  J»NP2)  +  .25*F  < J » NP1 ) 

75  CONTINUE 

C  CALCULATION  OF  AUXILIARY  PARAMETERS 

DO  80  I *2  » NP? 

TEMP < I )  =  (F(NH. I ) -U ( I )*U( I)/2.0)/CF  < I  ) 

IF  <TEMP( I ) .LE.O.O)  TEMP<I)=TINF 
TFMPE(  I)  *  TEMPI  I) 

AO  CONTINUE 

CALL  DENSTY 
DO  85  I  *  2 » NF'2 
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F0< J  '  =  R  H  0  < n *<53 . 3* TEMP < I ) +0 . 5*U ( I ) *U < I ) /32 . 2 ) / 1 4  4 . 
Hi.  AMACH  <  I  )=(/(  I  )  /SORT  (  1.4*53. 3*32 . 2*ABS  (  TEMP  (  I  )  )  ) 

90  YP  =  XU 

xu-xn 

ft  I=rPEI  +  PX*(R(  1  )*AMI-R<NF‘3)*AME) 

IF  < INTG. GT .  50  .AND.  NOFT  .EO.  1)  GO  TO  95 
GO  TO  105 

95  DO  100  I  =3,NP1 

IF  <ABS<DIFF<  I ) )  .GT.  1.5)  GO  TO  15 
00  CONTINUE 

THE  TERMINATION  CONDITION 
Of.  IF<XU.LT.XL)  GO  TO  25 
GO  TO  10 
10  STOP 
END 


SUBROUTINE  CONST 
REAL  MACH 

COMMON /BLK  2 /N  » NP 1 » NP2 » NP3 » NEQ » NF'H  .KEX.KIN.KASE.KRAD.PEI  »AMI  ,  AME  , 

1  DPDX.XU.XD.XP.XL.DX.  INTG  ,  CSALFA . HGS  .  C  INF  <  1  0  )  ,  PREF  <  10)  , 

2  PR< 10) »P< 10) ,U<050) »F< 10.50) . R < 050 ) . RHO < 050 ) » OM < 050 ) . Y<050) . 

3  DEN  <  050 ) » AMU  <  050 ) 

COMMON /BLK 3/ KRP . MACH.KR. INDI < 10) » INDE  < 10) . SUMR » AK . ALMG . BE T A . TAU I 

1  . TAUE.VINF.DELTA.TU. TINFfQS.YL.UMAX.UMIN.FR. YIP. YEM.TR.HU, 

2  HF< 10) , GAMA (10) » A J I < 10) »AJE< 10) , SC < 050 ) . AU < 050 ) . BU<050>  »CU<050 ) , 

3  A  < 10.050) ,B< 10.050) »C< 10.050) .TEMPE< 050) .TEMP (050) .P0< 050) . 

4  AMACH<050) »CP<050> 

COMMON /BLK 4 /NH . NNO . NNF .N0.NN.N02.NN2.NE  » I ND ( 1 0  )  . STO . AKS . RT . F  T . 

1  AMT  »  DRAGC . STNO . T AU I U , QU . UGU »  UGP 

C0MM0N/BLK5/KRI .PINF ,DIFF<50 ) . NOPT , NSK AT 
COMMON/CUALL/CU 

COMMON/SANJ/SIGMEX.SIGMAB.F' DEN, SIZE 
UNITS  *****  TEMPI  R  )  **•*  PINF  <  Ll«F/ IN2  )  ************** 

UNITS  *****  DEN <  L BM/F T 3 )  **  AMU < L BM/F T -SECO  **  CP < F T 2/SE C 2 - R > 

IF  RR* 1  INTEGRAL  RADIATION  ANALYSIS  IS  USED 
If  M *2  RADIATION  EFFECTS  ARE  NEGLECTFD 
IF  K  R*  3  OPTICALLY  THIN  ANALYSIS  IS  USED 

URITf (A.*) ' INPUT  THE  VALUE  T INF . TU. PINF .MACH. XL »KRI .NSKAT , 

1  CU.  (IS, REFT  ,SI  GME  X  ,  S I GMAB  .  8 1  ZE  » PDCN  ' 

MR1  Tf M  .*)' INPUT  THE  VALUE  TINF  »  T  W , P I  NT  .MACH, XL  »KF1 . WSK  AT • 

1  CU.OS.REF  T .SIGMEX.SIGHAB.SIZF .PDCN 


REAr*(5.*)TINF, TW, P INF, MACH. XL.KRI , NSKAT , CU , OS » REFT , S I GhE X , 

1  SIGMAB.SIZE  ,  PDEN 

UR I7E< 1 ,*>TINF,TW,PINF, MACH, XL.KRI .NSKAT.CW.QS.REFT.SIGMEX. 
t  SIGMAB.SIZE, PDEN 
KR  =  KRI 
TR-537.0 
X(J=  0.03 

DINF=PINF*144./< (53.3) *T INF) 

UISCM 1 .2267E-05)*(TINF/TR)**.76 
VINE  =  MACH*SGRT<1. 4*53. 3*32. 2*TINF) 

UGD  =  VINE 

DELTA=0.37*XU/< <VINF*DINF*XU/VISC)**.2> 

I N  [i  (  1  )  =  1 
NH  =  1 
AK= . 435 
ALMG  = . 09 
FR* . 0 1 
PREF  < 1 ) =0 . 9 
PREF ( 2 ) =0 . 9 
PREF <  3 ) - 0 . 9 
PR ( 1  )  *0 . 7 
PR ( 2  )  *0 . 7 
PR ( 3  )  =0 . 7 

P  < 1 ) =  3. 68* ( PR ( 1 ) /PREF  <1)-1.0)*((PR( 1 ) /PREF ( 1 ))**(-. 25) ) 
P(2)*5.0*(PR(2)/PREF(2)-1.0)*( ( PR ( 2 ) /PREF < 2 ) )**(-0.25> ) 
P<3)*5.0*<PR<3)/PREF<3)-1 .0)*( ( PR (3)/PREF(3) )**<-0.25) ) 
RETURN 
END 


SUBROUTINE  BEGIN 

COMMON /Bl K  2/N »  NP 1 f NP 2 . NP 3 » NEO t NPH r KE X » K I N r RASE  f KRAI* . PE  I t AMI » AME . 
1  PPDX.XU.XD.XP.XL  .PX.INTG.CSACEA.HGS.CINE  (  10)  »  PREF  <  10)  . 
i  PR<  10)  »P(  10)  f  U  <  050 )  t  F  (  10.050)  »R<050>  ,RM0<050)  ,0M<050>  .Y(050>  . 

3  DEN ( 050 ) »  AMU<  050 ) 

COMMON /BL  K  3/KRP »  MACH  *  KR . I NP I ( 1 0 ) . INPE <1 0 ) » SUMR , AK , ALMG , RE T A , 

1  TAUI »  TAUE .VINE .DEL  TA. TW. TINE  »  OS  »  YL »  UMAX  »  UM IN » ER  »  Y IP »  YEN  »  TR ♦ MW  » 

?  ME  < 10 ) .GAMA*  1 0 ) . AJ1 ( 1 0 ) . A JE  <  10  >  »  SC  < 050 ) . AU ( 050) »  RU<  050 ) . 

?  CU  *  050 > . A*  1 0 . 050 ) . B (  10.050 ) »  C ( 1 0. 050  > . T  E  MFT  < 050 ) . TEMP <  050 ) • 


(,<* 


4  PO(OSO) »AMACH(050> .CP(050) 

COMMON/BL K4/NH.NN0. NNP. NO , NN . N02 > NN2 . NE. IND( 10) r STO , AKS t RT , FT , 
J.  AMT.DRAGC.STNO.TAUIW.GW.UGU.UGD 

C0MM0N/BLN5/KRI .  PI NF  ,  D I FF  <  50  )  .  NOPT  ,  NSKAT 

COMMON/CWALL/CW 

DIMENSION  CINW<10) 

C  PROBLEM  SPECIFICATION 
KRAD  =  0 
KIN  ~  1 
NEQ  =  4 
KEX  -  2 
KASE~2 

IF ( KIN . EG  < 1 . OR . KEX . EQ . 1 )  KASE  =  1 

NPH=NEQ- 1 

NF‘1  =N+  1 

NP2=N+2 

NP3=N+3 

C  INITIAL  VELOCITY  PROFILE 

U< 1 )=0.0 
U( 2) =0 . 0 
DO  10  I  =3 » NF'l 

Y(  I )=DELTA*( (FLOAT! 1-2) )/FLOAT(N>  ) 

10  U< I)=VINF*(Y(I)/DELTA>**.143 
Y ( NP3 )  =  DELTA 
U ( NP3 ) =  VI  NF 

C  CALCULATION  OF  SLIP  VELOCITIES  AND  DISTANCES 

BETA* . 1 43 

GO  TO  (15.20.25)  .KIN 
15  U(2)=U(3)/(1.+2.*BETA) 

Y(2)=Y(3)*BETA/(2.+BETA) 

GO  TO  35 

20  U1 1*U( 1 >*u< 1 > 

U13*U< 1 )*U(3) 

U33*U ( 3 ) *U  ( 3 ) 

SQ=84.*UU-12.*U13  +  ?.*U33 

U<2)*< 16.*U11-4.*U13+U33)/(2.*(U( 1 >+U(3) )+SQRT(SG> ) 
Y(?)*Y<3)*<U<2)+U(3)-2.*U(1))*.5/(U(2>+U(3)+U(1)) 

GO  TO  35 

25  IF ( KRAD «  NE • 0  )  GO  TO  30 
U<  2>«(4 . *U( 1 )-U( 3 ) >/3. 

Y ( 2 ) *0 . 

GO  TO  35 
TO  U<  2>*U<  1  ) 

Y(2)*Y(3)/3. 

GO  TO  (40.45.50)  .KEX 
40  U<NP2)*U(NP1 )/( 1 ,+2.*BETA) 

Y(NP2)*V<NP3)-< Y(NP3)-Y(NP1 ) ) *BETA/ < 2 . +BETA > 

GO  TO  55 

45  U1 1 *U<NP1 )*U<NP1 ) 

U1 3*U<NF1 )*U<NP3) 

U33«U(NP3)*U<NP3> 
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SQ=84 .  CU33- 1 2 .  tUl 3+9 .  *U  1 1 

U  (  NF‘2 )  ~  (  16.  *U33- 4.  *U13+U1  1  ) / ( 2 . * < U < NP 1 ) +U < NP3 )  )+ SORT (SO)  ) 

Y  <  NP2  > =Y  < NP3 ) - ( Y  <  NP3 ) - Y ( Nf 1 ) ) * ( U ( NP2 ) +  U ( NP 1 ) -2 . *U ( NP3 ) ) * . 5/ 

1  ( U <  NP2 >  +U  ( NP 1 ) +  U  <  NP3 ) ) 

GO  TO  55 

50  U(NP?>=(4.*U<NP3>-U<NP1 ) )/3 

Y ( NP2 ) = Y  <  NP3 ) 

55  CONTINUE 

INITIAL  PROFILES  OF  OTHER  DEPENDENT  VARIABLES 
CP  1 =5990 . 0 
CP2=5990 . 0 
F(2»1)=1.0-CW 
F ( 3 » 1 ) =CW 
F ( 2 » NP3 )  =  1  .0 
F  <  3 »  NP3 ) =0 . 0 
TEMP  <  NP3 ) =  TINF 
HCU  =  0.0 
HCG  *  0.0 
HGS=CP1*TINF 

HW=  (CP1*F(2* 1 >+CP2*F<3. 1  )  >*TW 

HG=HGS+0.5*VINF*VINF 

F ( NH » 1 ) =HW 

F  ( NH »  NP3 )  =  HG 

DO  60  1  =  3 »  NP 1 

F ( NH » I ) =  HU+(HG-HW)*( Y( I > /DELTA)**. 143 
F<3.I)=CU*(1.0-U<I>/VINF) 

F  <  2 » I )  =  1 • 0-F ( 3 . 1 ) 

60  CONTINUE 

CALCULATION  OF  CORRESPONDING  SLIP  VALUES 
DO  100  J=  1  »  NPH 
GAMA (  J )  =  • 143 
GO  TO  (65.70.75) .KIN 

65  F  <  J »  2 ) =F  (  J » 1  )  +  (  F  <  J  »  3  )  -F  (  J  .  1  )  )  ♦  ( 1  • +BETA-GAMA ( J ) >/( 1 . +BETA+ 

1  GAMA  <  J  > ) 

GO  TO  80 

70  G=(U(2)+U(3)-8.*U(1 ) )/(5.*(U(2)+U(3> >+8.*U( 1 ) ) 

GF=(1.-PREF(J))/(1,+PREF( J) ) 

GF=(G+GF)/(1.+G#GF) 

F< J.2)=F( J,3)*GF+( l.-GF)*F( J.l  ) 

GO  TO  80 

75  F(  J.2)*F( J. 1  > 

IF(KRAD.EQ.0)F( J»2)=(4.*F( J»1 )-F( J.3))/3. 

80  GO  TO  (85.90.95) .KEX 

85  F< J»NP2)=F( J»NP3)  +  (F( J.NP1 ) -F ( J . NP3 ) ) * ( 1 . +BETA-GAMA ( J ) > /<  1  .  + 
1  BET A-GAMA  <  J ) > 

GO  TO  100 

90  G=(U(NP2)+U(NF1 )-8.*U(NP3) ) / ( 5 . # ( U ( NP2 ) +U ( NP 1 ) )+8.*U(NP3> > 

GF=  <  1  .  -F'REF  (  J  )  )  /  ( 1  ,  +PREF  <  J  )  ) 

GF=(G+GF>/( 1 ,+G*GF> 

F ( J  »  NP2 )=F ( J.NP1 )  #GF+ ( 1 . -GF ) *F ( J , NP3 ) 

GO  TO  100 

95  F< J»NP2)=(4.*F( J.NF3)-F( J.NP1 ) )/3. 

00  CONTINUE 
RETURN 


SUBROUTINE  READY 

COMMON /BLK2/N, NT  1 , NP2 » NP3 , NEO » NPH , KEX , N IN » KASE » KRAD » FE I  *  AM  I ,  AME, 
1  DPDX.XU.XD.XP.Xl  ,  DX »  INTO  ,  CS  AL  FA  ,  HGS  »  C I NF  (  10)  ,PREF(  10)  , 

?  PR (10) ,P< 10) ,U< 050) *  F ( 1 0  »  050 ) ,R(050>  ,  RHO ( 050) »0M(050) ,  Y( 050) , 

3  DEN (050) .AMU (050) 

COMMON/BLK  3/KRP  »  MACH  »  KR . INDI ( 10) . I N  D  E ( 10) . SUMR , AK . ALMG . BET A , TAU I 

1  .TAUE.UINF, DELTA, TU.T INF, GS.YL, UMAX, UM IN, FR, YIP, YEM,  TR.HU. 

2  HF(10) .GAMA ( 10).AJI<10)»AJE(10).SC<050) .AU(050) »BU(050) »CU(050) , 

3  A  (  10,050)  ,B(  10,050)  ,C(  10.050)  ,TF.hPE(050>  .TEMP  (050)  ,PO(050)  , 

4  AMACH(050) .CP(050) 

C0MM0N/BLN4/NH,NN0»NNP»N0»NN»N02»NN2»NE» IND(IO) .STO.AKS.RT ,FT  » 

1  AMT.DRAGC.STNO.TAUIW.QU.UGU.UGD 
CALL  RAD(XU»R( 1 ) .CSALFA) 

C  Y  NEAR  THE  I  BOUNDARY 

GO  TO  (10, 15,20)  .KIN 

10  Y ( 2 ) = ( 1 . +BETA)C0M(3)C4./( ( 3 . *RHO ( 2 ) +RHO ( 3 ) ) * ( U ( 2 ) +U ( 3 ) ) ) 

GO  TO  25 

15  Y(2)=12, *0M(3)/( ( 3 ♦ CRHO ( 2 ) +RHO  < 3 ) ) C ( U ( 2 > +U < 3 ) +4 . CU ( 1 >  >  > 

GO  TO  25 

20  Y(2)  =  .5C0M(3)/(RH0( 1 )CU( 1 >  ) 

25  Y(3)=Y(2M  .  25C0M(3  )C(1./(RH0(3>CU(3> ) +2 . / ( RHO ( 3 ) CU ( 3 ) +  RHO ( 2 ) * 

1  U  ( 2 )  )  ) 

C  Y  S  FOR  INTERMEDIATE  GRID  POINTS 
DO  30  1  =  4  ,  NF'l 

30  Y( I >=Y( 1-1 )+.5*(0M( I )-0M( 1-1 ) )*( 1 . / ( RHO ( I ) CU ( I ) ) + 1 . / ( RHO ( I - 1 ) * 

1  U(I-l))) 

C  Y  NEAR  THE  E  BOUNDARY 

Y(NP2)=Y(NP1 )+.25C(0M(NP2)-0M(NPl ) ) * ( 1 . / ( RHO ( NP 1 ) *U ( NP 1 ) )+2./ 

1  (RH0(NP1 )CU(NP1)RRH0(NP2)CU(NP2) ) ) 

GO  TO  (35,40,45) ,KEX 

35  Y(NP3)=Y(NP?)+( 1 . +BETA ) C (0M(NP2 > -0M( NPl ) )C4. /( ( RHO ( NP 1 ) +3 . * 

1  RHO  (  NP2  )  )  *  (  U  ( NF‘l  )  +U  (  NF‘2  )  )  ) 

GO  TO  50 

40  Y<NP3)=Y(NP2)+12.C(0M(NP2)-0M(NP1 ) )/( ( RHO ( NP 1 ) +3 . CRHO ( NP2 ) >C( 

1  U(NP2)+U(NP1 )+4.*U(NP3) ) ) 

GO  TO  50 

45  Y(NP3)=Y(NP2)T.5*(0M(NP2)-0M(NP1 ) ) / ( RHO ( NP3 ) *U ( NP3 ) ) 

50  IF (CSALFA. EO.O. .OR.KRAD.EG.O)  GO  TO  60 

DO  55  1=2, NP3 

55  Y( 1 >=2.CY( I ) CPE I / ( R ( 1 ) ♦  SORT ( R ( 1 ) CR (  1 > +2 . C Y ( I ) CPE  I CCSALFA ) ) 

GO  TO  70 

60  DG  65  1=2. NP3 
65  Y( I )=PEICY( I )/R( 1 ) 

70  Y(2>=2.*Y(2)-Y(3> 

Y(NP2)-2.CY(NP2)-Y(NP1 ) 

C  CALCULATION  OF  THE  RADII 

DO  75  1=2, NP3 
IT (RRAD.EQ.O)  R(I)=R(1) 

IF  (KRAI'.NE  .0  )  R(  I  )=R(  1  >+Y(  I  )  CCS  ALFA 
75  CONTINUE 
RETURN 
END 


SUBROUTINE  DENST Y 

C0MM0N/BLK2/N,NP1  »NP2»NP3,NEG,NPH,KEX,KIN,  KASE , KRAD , FE I  ,  AMI  ,  AME, 

1  DPDX.XU.XD,  XP,  XL.DX,  INTG , CSALFA , HGS » C INF  (  10)  , PREF  (  10). 

2  F'R  (10)  .  P  ( 10)  ,  U  <  050  )  ,  F  (  10.050)  »R(050)  »RH0(050)  ,0M(050)  , Y(050)  , 

3  DEN(050) ,AMU(050) 

COMMON/ BLK3/KRP, MACH, KR, INDI (10) , INDE( 10) . SUMR , AK » ALMG . BETA » TAU I 

1  f TAUE.O INF .DELTA. TU.TINFf OS. YL. UMAX. UMIN.FR, YIP, YEM.  TR.HU, 

2  HF( 10) ,GAMA( 10) , AJI ( 10) , AJE( 10) , SC ( 050 ) » AU ( 050 ) , DU (050) ,CU(050 ) , 

3  A ( 10,050) ,B( 10,050) , C < 1 0 , 050 ) , TEMPE ( 050 ) , TEMP ( 050 ) , PO ( 050 ) . 

4  AMACH(050) ,CP(050) 

C0MM0N/BLK4/NH»NN0»NNP,N0,NN»N02,NN2»NE, IND( 10) » STO , AKS , RT , FT , 

1  AMT.DRAGC.STNO.TAUIW.QW.UGU.UGD 
DO  115  1  =  1 »  N P 3 

RHO ( I ) =  DEN< I )*<TR/TEMP< I ) ) 

115  CONTINUE 
RETURN 
END 


SUBROUTINE  RAD ( X . R 1 » CS ALFA ) 

APPLICABLE  TO  PLANE  FLOWS 
CSALFA  =  1.0 
R 1  =  1.0 
RETURN 
END 

FUNCTION  UISCO(I) 

C0MM0N/BLK2/N.NP1 , NP2 , NP3 , NEG , NPH , KEX , K I N , KASE » KRAD , PE  I , AMI , AME. 

1  DPDX,XU,XD,XP,XL,DX,INTG,CSALFA,HGS»CINF( 10) ,PREF< 10) , 

2  PR< 10) »P( 10) »  U ( 050 ) ,F< 10,050) »R(050) , RHO (050 ) ,0M(050) , Y( 050 ) , 

3  DEN(050) ,AMU<050) 

C0MM0N/BLK3/KRP, MACH, KR, INDI ( 10) »INDE( 10) .SUMR, AK, ALMG, BETA, TAU I 

1  ,TAUE,UINF,DELTA,TW,TINF,QS»YL,UMAX,UMIN,FR, YIP, YFM,  TR.HU, 

2  HF( 10) ,GAMA( 10) ,AJI ( 10 ) , AJE( 10) , SC (050) , AU(050) , BU(050 ) ,CU(050 ) , 

3  A  (  10,050) , B( 10,050) »C( 10,050) , TEMPE (050) , TEMP (050) ,P0(050> . 

4  AMACH(OSO) . CP (050) 

C0MM0N/BLK4/NH»NN0»NNF‘,N0»NN»N02»NN2,NE»  IND(  10)  ,STO,AKS,RT,FT. 

1  AMT , DRAGC.STNO, TAU I W , QU , UGU , UGD 
VISCO=AMU( I )*(TEMP( I )/TR)**. 76 
RETURN 
END 


SUBROUTINE  I.  ENG  T 

COMMON/HLKr’/N.  NFl  ,NF'2.Nr3.Nf  Q.NRH.KEX.KIN.KASE  .KRAD.REI  .AMI  ,AME  . 
1  [irnX.XIJ.Xh.Xf.XL.nX.INTG.CSA!  FA.HGS.CINFf  1  0  >  .  FREE  <  1  0  )  • 

;>  PR(  ?  0>  »!•'  1  0  )  .  U(  050  )  .  F  (  )  0  .  Af,0  )  .  R  <  050  >  .  RHO<  050  )  .  OH  (  050  )  .  Y  <  050  '  , 

3  OF N ( 050 ) »  AHU ( 050 ) 

COMMON/ Bl  K  3/KRP. MACH.KR , INDI < 1 0 ) . INDE  < 1 0  > . SUMR. AK . AL  MG  »  BETA . T AU ! 
1  »  T  AUE  »  0  1  NF  ,  Itfl.  TA.TU.T  I  NF  •  PS  .  Y  L  ,  UM  AX  ,  UH I  N  ,  FR  ,  Y  1  F  ,  YE  M  ,  TR  .  HU  . 

?  HE  < 1 0  > . GAMA (  1 0  > »  AJI < 1 0 ) • A JF <  1 0  > • SC ( 050 ).AU( 050 ) »  BU ( 050 ) . CU ( 050 '  . 

3  A( 10 . 050 ) . B< 10.050) »  C<  1 0.050) »  TEMF  E  <  050 ) . TEMF  <  050  > ,F0<  050 ) . 

4  AMACH<050) »CP<050) 

COMMON/KL  K4 /NH .  NNO  .  NNF'  ,N0*NN.N02.NN?.NE.  I  ND  <  10)  .STO.AKS.RT, FT  . 

1  AMT.DRAGC.STNO.TAUIU.QU.IJGU.UGD 
C  Sf  ARCH  rOR  MAX  T  HUH  AND  MINI MUH  VEIOCITITES 
10  UMAX  =  U ( 1  ) 

U  M  I  N  -  U  (  1  ) 

DO  15  J- 3  *  NF  3 

IE( J.EQ.NP2)  GO  TO  15 

IF<U' J) .GT .UMAX)  UMAX  =  U ( J  > 

IF<U< J) .LT.UMIN)  UM I N  =  U  <  J  > 

15  CONTINUE 

DIF=ABS(UMAX-UMIN ) *FR 
C  SEARCH  FOR  THE  I  BOUNDARY 

IFfKIN.NE.2)  GO  TO  35 
U? 1 = ABS  <  .5*(U<?>+U(3) )-U< 1 >  ) 

IF(U21 .LT.PIF)  GO  TO  20 
YIP=SORT( DIF/U21 )*<Y(2)+Y(3>>*.5 
GO  TO  40 
20  Jr  2 
25  J- J41 

u  J  i  =  u  <  n  -  u  ( i ) 

IF(ABS(UJ1 ) .GE.PIF)  GO  TO  30 
GO  TO  25 
30  A 1  =  1  . 

IF ( UJ1  . L  T . 0  .  )  Alr-i. 

YIR*Y<J-1)*<Y<J)-Y<J-1>)*<U<1>4A1*DIF-U(J-1>>/<U<J>-U<J-1>' 

GO  TO  40 
35  YIF=0. 

C  SEARCH  NEAR  THE  E  BOUNDARY 
40  IF ( K E X . NE . 2  >  GO  TO  60 

!J2 1  *  ABS  <  . 5t <U<NR1 >4U< NF 2 ) > -U< NF  3  >  » 

IF(U21 .LT.DIF)  GO  TO  45 

YE M* SORT ( DIF/U21 )*  <  .5* ( Y ( NF 1  ) ♦ Y  <  NF  2  >  > - Y ( NR3 )  >  4  Y ( NF  3 ' 

GO  TO  45 
4  5  J  =  NF  2 

50  JrJ-1 

UJ1 =  U<  J  >  -IM  NF  3  > 

IF  (ABS(U.)l  >  .G£  .DIF  >  GO  TO  55 
GO  TO  50 
55  Ml  r  1  . 

I  F  <  U  1 1  .  l  T  .  0  .  )  A  1  -  -  1  . 

YF  M-  Y  (  J4  1  *  4  (  Y  l  .1  >  -  Y  <  I  4  1  '  )  #  (  U  ‘  NF  3  »  4  A  1  •  D  I  F  I  I  ■  l  4  .  '  i  <  U  1  J  '  U  '  1  *  1 

GO  TO  65 
4.0  y  F  M-  V  (  NF 1  » 

4.5  Yl  *Y[H  Y  IF 
Rf TURN 
ENp 


’4 


SUBROUTINE  MASS  ( XU »  XP • AM ) 
r  AFT LTCABLE  TO  F'OROUS  UAL  l 

AM  -  0.0 
RETURN 
END 


SUBROUTINE  ENTRN 

COMMON/BLF.  2/N»NF’l  ,  NP2  ,  NF'3  ,  NEQ  ,  NPH  ,  KE  X  ,  K I N  .  K  ASE  ,  KR  AD  ,  PE  I ,  Ah  1 ,  AME  , 

1  PPPX  ,XU»  XP.XP.  XL  ,I»X.  INTG.CSALFA.HGS»CINF(  10)  .  PREF  (10)  » 

2  PR  (  10) ,F( 10) ,U(050)  ,  F  ( 1  0  ,  050  )  ,  R  (  050 )  , RHO ( 050 >  » OM ( 050 )  ,Y(050) , 

3  DEN ( 050 ) »  AMU ( 050 ) 

COMMON/BLK3/KRF'»MACH»KR»INr*I(10)  ,  INDE( 10)  ,  SUMR  ,  AK  »  ALMG  ,  BETA  ,  TAUI 
1  .TAUE, VINE, DELTA, TU , T I NF , GS » YL . UMAX » UMI N , FR , Y I P » YEM ,  TR.HW, 

?  HF( 10) »  GAMA ( 10) ?AJI ( 10) »AJE( 10) » SC (050) ,AU(050) , BU(050) »CU(050) , 

3  A  (  10 r 050) .»( 10.  050 )  » C  ( 1 0 , 050 ) » TEMFE ( 050 ) r TEMP (050) »P0(050) , 

4  AMACH(050) .CFK050) 

COMMON/BLK4/NH,NNO,NNP,NO.NN,N02,NN2,NE. I  NEK  10) , STO , AKS » RT r FT , 

1  AMT,PRAGC»STNO»TAUIU»GU,UGU»UGP 
C  THi 5  SUBROUTINE  USES  THE  M I  X I NG-LENGT  HYPOTHESIS. 

GO  TO  (70,75.80) .KIN 
70  GO  TO  85 

75  AM  I =8. ♦  R  H  0  ( 1 >*( (ALMG*YL)/(Y(2)+Y(3) ) ) **2*ABS ( U ( 2 ) +U ( 3 ) -2 . *U ( 1 ) ) 

GO  TO  85 
80  AM  1=0. 

85  GO  TO  < 90.95. 100) ,KEX 
90  RETURN 

95  AME  =  -8.*RH0(NF3)*( ( ALMG* YL > / ( Y ( NP 1 ) + Y  (  NP2 ) -2 . *Y ( NP3 ) ) )**2*ABS( 

1  U(NP1 )+U(NP2)-2.*U(NP3)  ) 

RETURN 
100  AME  =  0 . 

RETURN 

ENB 


SUBROUTINE  WALL 

COMMON /BLK2/N.NP1 , NP 2 , NP3 , NEQ , NPH , KE X , K I N » K ASE , KR AD » PE I » AM  I » AME • 
1  PF  DX . XU. XP, XP, XL . PX, INTG » CSALFA , HGS  *  CINF ( 1 0 ) »  PREF ( 10), 

?  P  R  <  10*  »  F'  <  10)  ,U<  050)  .  F  (  10,050)  ,R(  050)  ,RHO<  050  )  ,0M<  050)  ,  Y  (050)  , 


1 


COMMON/ PL K3/KRP, MACH *KR, I NP I ( 1 0 ) . I N PE < 1 0 > » SUMR * AN , AL MG . PE T A , TAUI 

1  .  TAUE*  U  INF,  DELTA,  TU »  T  I NF  »  OS »  Yl.  *  UMAX.UM  IN *FR, YIP,REM,  TP  .  HU  • 

2  HF( ] 0> *  GAMA ( 1 0  >  » A JI < 1 0  >  *  A JC <  J  0 > » SC (050) * AU<050) ,  PU<  050 ) *CU(050) . 

3  A  < 10*  050 ) *B<  10*050) »C<  10*050) »TEMFE(050) . TEMP( 050 ) »P0( 050 ) . 

A  AMACH ( 050 ) *  CP ( 050 ) 

COMMON /PI.  M/NH  *  NNO.NNP.NO  ,  NN  »  N02  *  NN2  *  NE  *  IND < 1 0 >  » STO . APS  * RT , F T , 

1  AMT*nRAGC*STNO*TAUIU*GU,UGU*UGD 
C  CAl.CUl.  AT  I  ON  OF  BETA  FOR  THE  I  BOUNDARY 
I F ( KEX . NE . 1 )  GO  TO  25 
10  YI=Y(NP3>-.5*( Y<NP1 )+Y(NP2)  ) 

UI=.5*(U(NP2)+U(NP1 ) ) 

RH®,25*(3.*RH0(NP2)+RH0(NP1 ) ) 

RE=RH*UI*YI/UISC0<NP3> 

FP-DPBX*YI/(RH*IJI*UI ) 

AM=AME/(RH*UI ) 

CALL  UF 1  (  RE  >  FF‘  *  AM  *  S ) 

BET A® SORT <ARS(S+FP-fAM) )/AK 
TAUE=S*RH*UI*UI 
IF(NEG.EG.l)  GO  TO  20 

C  CALCULATION  OF  GAMA  S  FOR  THE  E  BOUNDARY 
DO  15  J=1»NPH 

CALL  WF2  <  RE  *  FP  *  AM  *  PR ( J ) *  F‘REF  (  J  )  *  P  (  J )  *  SF  ) 

GAMA< J)=(SF+AM)*PREF( J ) / < AK*AK*BETA > 

IF(INDE( J) .EG. 1 )AJE< J  )  =SF*RH*UI * ( F ( J*NP2)+F( J,NP1 ) -2 . *F ( J * NP3 ) )* 
1  .5 

15  CONTINUE 
20  IF ( K I N . NE . 1  ) RETURN 
C  CALCULATION  OF  PETA  FOR  THE  I  BOUNDARY 
25  YI=.5*<Y(2)FY<3> ) 

UI  =  .5*<U(2)+U<3>  ) 

RH®,25*(3.*RH0(2)+RH0(3) ) 

RE=RH*UI*YI/UISCO< 1 ) 

FP=DPDX*YI/(RH*UI*UI ) 

AM  =  AMI/(RH*UI  ) 

CALL  UF 1 ( RE  *  FP  » AM  *  S ) 

PETA-SQRT(ABS(S+FP+AM) )/AK 

TAUI=S*RH*UI*UI 

TAUIW  =  TAUI/32.2 

DRAGC  =  (TAUI/(RH0<NP3)*U<NP3)*U<NP3) > >*2. 

IF(NEG.EO.l)  RETURN 

r  CALCULATION  OF  GAMA  S  FOR  THE  I  BOUNDARY 
DO  35  J= 1  *  NPH 

30  CALL  UF2  <  RE  *  FP  *  AM  *  PR  <  J ) • PREF ( J ) »  P ( J ) »  SF ) 

HI®  .5»<F(NH* 2)+F(NH*3)  ) 

PHI  *  (UI*UI/(2.0*(HI-F<NHf  2)  )  )  )*(  (  1  .  O-F'REF  (  NH  )  )  *SF*PREF  (  NH  )  / 

1  <  AK*AK*BETA)-2.0#BETA#( 1 . O-PREF (NH) ) ) 

IF( J.EO. 1 )  STN0=SF*(RH*UI/<RH0(NP3)*U(NP3> ) ) 

GAMA( J) ®(SFFAM)*PREF( J ) / ( AK*AK*BETA ) -PH  I 

IF< INDI ( J) .EG. 1 ) AJI< J)=SF*RH*UI*(2.*F( J* 1 )-F( J* 2)-F( J*3) )*  .5 
•5  CONTINUE 

T*0  r  0.015 

)*>  -  TK0*(0. 5t  (  TEMF’(  3  ) +TU  ) /TR  )  **0 . 76 
»W  *  T**<  TEMP' 3>-TU)/(  Y(3)*3600.  ) 
fc  r  Tiji.N 
f  **  f ' 


*  ■ 


’  I  N(  y  I  '  .«  .  aa  ,  < 

'  -  H  ►  .I'KI  .  *1A(H.K  .  ]  W  l<  !  •  J  0  •  .  :  *•  1  ♦  •  I'1  .  S'  'AF  .  A,  .A)  At.  .  M  U.  .  'A 

«  T  A  1  It  ,  V  1  *  I  .Ml  1  f.  .  T  u  .  T  I  Ml  ,  f|  *  .  *  I  .  1 1 A  A  »  .  I A  !  A  .  ♦  I-  .  »  1  »  .T|A.  »  »■  ,  M  U  . 

Hf  «  t  0  »  .  G  AHA  ‘  1  <)  >  .  A  ;  !  *  1  0  •  .  A  l|  1  1  A  *  .  M  -  O*  >j  •  .  At)  <  0*,0  .  f* I  •  •  O  V*  .in.  ■  . 

A  «  1  A  .  'I  *.  A  •  ,  f.  l  1  n  ,  .  «  0  .11’.  .  O'.  .  '  t  «l  I  ■  •  *  <  .  *  »  f»  I  ■  '  •  I  I  .  I  '<  *  f 

AHAf  M  -  <  ',0  .  '  I  '•*,*) 

MJHHOH  /  Rf  M/MH.NNiJ.NNI  .NO.NN.NIJ.'.NN,  .  Nf  .  1  M||  .  1  <  !  1 1  ,  fl  |  >  . 

A  h  T  .  I •  ►  aG<"  .  c,  t  yn  .  r  ah  l  w ,  uy  ,  ik.ii  .  hg  [ 
af S  ~  AF 
RT  ^R*AFr 

‘jT  =  l  .  /  k  T  -  .  | '(  1  •  F  T  t  •  >  .  4  * •  ♦  .  0  6  7  :  <  •  F  *  •  *  •  1  >  ♦  .  0  J  /  1  !  •  F  ■  •  •  !  I- 

sto-  st 

IF  (  F  .Ell  .  0.  >  GO  TO  1 
F  t -F /  Af  0 

F  M  T  1  .  -  4  .  •  F  T  •  R  T  >585. ♦FT**?.  *,>•*. 4 
I F  <  F  M . I  T.0.1FH-0. 

ST=ST*F«*»J .A 
GO  TO  20 

IF(AM.FQ.O.  »  GO  TO  ?0 

AMT*AH/AFS 

AMM~  1  .  -  AMT  /  <  7 . 744R  T»*  <  -  1  .  1  7  >  ♦  .  956**>  T  »•  <  .  25  > 

ST=ST»AMM*»4 

S=ST*AKS 

RETURN 

END 


SUBROUTINE  UF2  <  R  *  F  » AM ,  PR . PRT , F . S ) 

C0MM0N/BLK3/KRF, MACH»KR» INDI< 10) . INDE (10) »  SUMR  ,  AK , AL MG . RE T A , T AU I 

,TAUE. VINF.DELTA, TU.TINF.GS. YL.UMAX .UMIN.FR. Y IP. YEM. TR.HU, 

HF( 10) ,GAMA< 10 ) , AJI < 10 ) , AJE( 10) , SC  <050 ) . AU( 050 ) , BU( 050 ) . CU<  050 )  . 
A  <  10,050)  .  B(  10,050)  »C<  10,050)  ,  TEMF'E  (  050  )  ,  TEMP  (050  )  ,PO(050  )  . 
AhACH(OSO) ,CP(050) 

COMMON/BLK4/NH,NNO,NNP,NO.NN,N02»N2»NE , I ND ( 10) »ST0»AFS»RT .FT , 

AMT»DRAGC,STNO,TAUIW»GU,UGU»UGD 

ST1=ST0/(1.0+P*SGRT(ST0) ) 

IF(F.EQ«0.0 ) GO  TO  10 

SSEP=1 ,725#RT**<-.  3333)*  <F-F6.0>**< -1  .  165) 

FD-.25*FT*RT/< 1 . F.0625*RT) 

ST1=ST1*< 1  .0-FD)+FD*SSEP 

ST=ST1*PRT 

S=ST*AKS 

RETURN 

END 


it  ■  '  I  It’ 


i  t-  ii '  :  «•* 

1  i*.|  ^  A  i  .. 

'*»*"*  h  »  N.  Ml  |  .  Ml  .Ml  l.nl  II.MlM.M  1,1  |N.MM  ,  I  RAF . I  |  I  >  AH  I  •  AMI  , 

’•  I  >.»l  .il.ll  .  »  .  I<  »  .  |  M  f  (,  .  I  SAt  I  A  .  MC,3  .  I  I  Ml  <  1  O  >  .  I  Rf  I  i  10  1. 

•  *  •  .  i  I  o  ,  ,i  <  o\r.  .  . »  ■  |  o  .  nw,o  •  .  i.  >  030  )  .  i«n «  030  »  .  om  030  » .  r  <  030  >  , 

•  !-|  M<  ■>*,■')  .  AM  1 1  ■  03  'I  ■ 

OMMOM  HH  ►  I  .  MAI  m  .  -  I.  .  I  Ml*  1  -  I  o  I  .  I  M(|(  <  I  0  1  .  S  UMR  ,  AR  .  At  MG  .  M  T  A  .  T  A I !  I 

!  .IftlK  .  V|MI  ,(i(l  T  A  .  ru.  T  (Mf  .0^  .  ri  .  IIMA  X  .  UK  I  M  .  I  k  ,  V  (  F ■  .  r  f  M  ,  T  K  .  HU  . 

Ml  ■  I  0  <  .  f.  AMA  ■  I  0  )  .  A  I  I  <  1  0  I  .  A  )(  (  1  0  )  .  SC  <  030  >  .  AU  <  030  *  .  RU  (  030  )  .  CU  <  030  )  , 

1  A  <  1  0  .  '-30  >  .  I< '  I  '  .  030  >  .  C  '  1  0 . 030  >  ,  Tf  Mf  f  f  030  >  .  TF  Ml  <  030  t  .  PO  *  030  >  . 

4  AMA<  M  1  03c  .I  f  <  030  I 

I  MMMOM/H  R  4/MM.MMO.MMl  . NO • NM • NO: . NM? , N(  .  (HIM  (0)  .ST0.AR5.RT  ,  F  T  , 

AM  T  ,  fiA  A  r,  (  .  S  ’  MO  .  T  AU  I  M  ,  CJW  .  UGIJ  .  UGD 

I  OMMOM  /  f«  l  ►  3  R  I  .  I  I  Ml  t  M  I  r  <  030  >  t  NOt‘  T  .  NSK  AT 

o  '  Ml  M3  IOM  t  In  -  030  .  T  [It  I  030  > 

I I  '  t  M  »  G  .  Ml  .  1  >  GO  TO  20 

WA  I  T  I  1  t,  ,  !  3  *  T  I  Ml  .  T  W  .  MAC  H  .  RRAl»  iKlN.NEOilfl.N.  MGS  »  My  »  PR (  J  )  » 

1  Vj  Ml  .  [.fi  TA.r  IMT  tMi  1 

WR  I  Tf  •  1  ,  1  3  >  T  IMF  .  TW  .MACH,  RRAD.MN.NEO.  REX  .N  .  MGS  »HU.  I'M  1  )  , 

(  V JNT  . DEl  T A ,P INF ,RR  I 

1*  » PfcMAT < /// 1 K , 7MT INF  »  FA.1/1X.5MTW  -  F6.1/IX. 

1  10MMACM  MO  *  F  4 .  1  /  l  X  »  7HKR AD  «  I2/1X.AMRIN  -  I2/1X 

.AMMtO  r  I  2 / 1 X • 6HRE X  -  I2/1X.4MN  -  I  3 / 1 X  *  AHHGS  *  E13.4/1 

'  X.3HMU  =•  I  1 3 . 4/ ] X . 1 3HFRANDTI  NO  *  F  4  .  1  /  1  X  *  7HV  I  NF  -  E13.4/1X. 

«  8H[>E  L  T  A  t  E13.4./IX. 'PINF  «  '.E13.4./JX.RR  *  'ID 

.’0  CONTINUE 

'  IF  INF .ME . D  GO  TO  25 

!F<F10AT(INTG-1)/30..NE.FL0MT((INTG-1)/50)>  RETURN 
23  CONTINUE 

Vise- ( 1 . 226 7E- 05)4 ( T I NF / TR >»*. 76 
DINF -PINFA1 44 . /! 33. 3*T IMF  > 

R  E  N  X  -  V  INF  4DINF  9XU/VISC 
DO  30  1-2.  HP? 

▼  I'L  <  i  )  -  r  <  I  >  /  r  <  ne  3  > 

30  TDl !  I  > -TEMPI  I >/TlNF 
I  F  I  R  R  .  E  Q  .  2  >  SUMR-0.0 
C  IF(KR.FQ.l)  CALL  RADWAl 

WRITE  < A. 35>XU.RENX»AHE»SUMR 
WRITE  < 1 . 33>XU.RENX.AME.SUMR 

T3  FORMAT ' 7M1  XU  *  2REI1.2.4H  FT.I0X.9H  RE . X  =  1PE9.2. 

I  4  X  * '  AME  =  ' .2PE13.4.5X. '  OR  *  ' . 2PE 13.4.  ' RTU/I SO-FT  SEC)  ') 

WRI Tf 16.40) DRAGC  » STNO .TAUIWrQW 
WRI TF (  l  ,  40  >  DRAGC .STNO.TAUIW.QW 

FORMATllOH  DRAGC  =  2PE  1  3 . 4 . 1  OH  STNO  =  2PE13.4.UH  TAUIW  = 
1  2FE1 3.4.9HLRF/S0-FT  » 0H  QW  =  2FE 1 3 . 4 . 1 5HDTU/ ( SG-FT  SEC)) 
OT--SUMRFQW 
QROW=-RUMR/QU 
WRITE! 6.4) 'THE  VALUE' 

WRITE! 1  ,  »  )  ' THE  VALUE  ' 

URITCC6.59>QT»GR0U 
WRITE! 1 »59>QT»0RQW 

59  FORMAT ! 7H  QT  -  2FE 1 3 . 4 . 1 5HDTU/ ! SO-F T , SEC >  . 

1  10H  QRQU  =  2PE13.4) 

IF!NOPT.EQ. 1 )  GO  TO  43 
GO  TO  75 
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DU  F  <  1  >M  (  If  MP!  141  >-TEMP!  I  >  )/TEMP!  I  )  >*100.0 
WRITE!6.55> 

UK  I  TE  <  1 .55) 

f  ORMAT  ! 1H0.8X . 1HY . 1 IX . 1HU. 1 1 X . 1HH. 1 1 X. 1HT  » 1  1  X  *  2HP0  *  10X. 1HM. 

10X . 3HRH0  »  8  X  * 4HDIFF / ) 

WRITE (6. 60) 

WRITE< 1 »60> 

FORMAT  < 6X. 4H< FT ) .7X.BH! FT/SEC ) .  4X» 1 OH < FT2/SEC2 ) .  4X  »  7H( DEG-R) .  4X 
9H!LRF/IN2>  » 14X.9H!LDM/FT3> .6X.9H! PERCENT)/) 

UR  I TE<  6.65) Y <NP3) .U! NP3 ) »F! 1 .NP3) 

WRITE ( 1 . 65  5  Y(NP3>  »U<NP3>  *F( 1 .NP3 ) 

FORMAT  < 1H  1P8E12.3) 

DO  70  J 1 =2  » NP2 
J?«NP2-J1E2 

WRITE ( 6.65) Y<  J2) »U<J2) »F< 1 »J2) »  TEMP ( J2 ) » 

P0<  J2 ) » AMACH ( J2 ) »RHO( J2) » F ( 2 » J2 ) » F ( 3 » J2 ) 

WRITE! 1 .65) Y(J2) »U<J2) »F< 1 »J2) .TEMPI J2) » 

P0< J2) .AMACH! J2) . RHO ( J2 ) . F ( 2 . J2 ) »F!3. J2) 

CONTINUE 

WRITE! 6. 65) Y! I ) »U! 1 ) »F! 1 . 1 ) 

WRITE! 1 .65) Y! I > ,U! 1 ) »F! 1 . 1 ) 

GO  TO  100 
WRITE! 6.80) 

WRITE! 1.80) 

FORMAT ! 1H0.6X. 1HY. 11X. 1HU .11X.1HH. 1 IX. 1HT. 1 IX . 2HP0 . 10X , 1HH , 

1  OX  »  3HRH0 . 10X.7HY/DELTA.10X.6HT/TINF/) 

WRITE <6. 85) 

WRITE! 1.85) 

FORMAT (6X.4H! FT) .7X.8H! FT/SEC ) .4X» 1  OH! FT2/SEC2) .4X.7H ! DEG-R ) »4X 
9HILDF/IN2) . 14X.9HILBM/FT3)/) 

WRITE! 6.90) Y!NP3> .U!NP3) .F! 1 .NP3) 

WRITE! 1 .90) Y!NP3) »U!NP3> »F! 1.NP3) 

FORMAT ! 1H  1P9E12.3) 

DO  95  J1=?,NP2 
J?=NP2- J1 +2 

URITE!6.90>Y! J2) »U! J2) »F! 1 . J2) .TEMP! J2) 

PO! J2) .AMACH! J2) » RHO! J2> .YDL! J2) .TDL! J2) 

WRITE! 1 .90) Y! J2) .U( J2) »F! 1 » J2) .TEMP! J2) 

PO! J2) .AMACH! J2) .RHO! J2> » YDL ! J2 ) .TDL! J2> 

CONTINUE 

WRITE! 6. 90) Y! 1 ) »U! 1  )  .F! 1 » 1 ) 

WRITE! 1 .90) Y! 1 ) .U! 1 ) .F! 1 . 1 ) 

CONTINUE 

IF(KR.EO.l)  WRITE! 6. 105) 

IF ! KR . EG . 1  )  WRITE! 1 , 105) 

F0RMAT!6X. 'RADIATION  EFFECTS  INCLUDED  -  INTEGRAL') 

IF ( KR . EG • 2 )  WRITE ! 6 . 1 10 ) 

IF! KR . EG . 2 )  WRITE(l.llO) 

FORMAT !6X. 12HN0  RADIATION) 

IF! KR • EG . 3 )  WRITE! 6. 115) 

IF!KR.EQ.3>  WRITE!1.115) 

F0RMAT!6X. 'RADIATION  EFFECTS  INCLUDED  -OPTICALLY  THIN') 

RETURN 


SUBROUTINE  COEFF 

COMMON/ BL K  2/N  *  NP 1  »  NP2  .  NP3  »  NEG »  NPH  . KEX . R IN .  RASE  »  KRAP  .  PE  I  .AMI  ,  AME  . 

1  DPDX.XU.XD.XP .XL .DX. I NTG . CSALFA . HGS . C I NF < 10) .PREF ( 10) . 

2  PR< 10) .P< 10) ,U<050> .F(  10.050) .R<050) .RH0<050) »0M( 050) . Y<  050 ) . 

7  BEN (050) .AMU <050) 

COMMON /BLK3/KRP .MACH.KR. IND1 ( 10) . I  NOE ( 10  > . SUMR . AK . ALMG . BE T A , T AU I 

1  .  T  AUE  .VINF.DELTA.TW.TINF.QS.YL.UMAX.UMIN.FR.YIP.YEM.TR.hu. 

2  HF< 10) . GAMA  < 10).AJI<10)»AJE(10).SC(050) »AU(050) . BU ( 050 ) . CU  (  050 )  . 

3  A ( 10.050 ) »B( 10.050) »C<  10.050) .TEMPE(050) . TEMP (050 ) »  P0( 050 ) . 

4  AMACH ( 050 ) » CP  <  050 ) 

COMM ON /BLK4/NH. NNO.NNP.NO. NN .N02.NN2.NE. I ND ( 1 0 > . STO . AKS . RT , F T , 

1  AMT.DRAGC.STNO.TAUIU.QU.UGU.UGD 

DIMENSION  G 1  ( 200 ) »G2( 200 ) . G3 ( 200 ) . D<  10.200 ) ► SI < 200 ) »S2  < 200 ) . 

1  S3  <  200 ) 

IF(KR.EQ.l)  CALL  RADPRM 
C  CALCULATION  OF  SMALL  C  S 

DO  10  1=2. NP1 
RA=.5#'R(I+1)+R(I)  ) 

RH= . 5*  v  RHO ( I  +  1 >  ERHO< I > ) 

UM=.5#<Umi)+U<I>) 

CALL  VEFF< I . I+l.EMU) 

10  SC< I )=RA*RA*RH*UM*EMU/(PEI*PEI ) 

C  THE  CONCESSION  TERM 

SA  =  R< 1 >*AMI/PEI 
SB=(R<NP3)*AME-R(1 >*AMI)/PEI 
DX=XD-XU 
DO  25  1=3. NP1 
OMD=OM< 1+1 )-0M( 1-1) 

P2= . 25/DX 
P3=P2/0MD 

P1=(0M( 1+1 ) -OM ( I ) )*P3 

P3=<  0H< I )-0M< 1-1 ) )*P3 

P2=3 . #P2 

Q-=S  A/OMD 

R2=-SB* . 25 

R3=R2/0MD 

Pi *- <0M ( 1+1 ) +3 ♦ *0M( I ) ) *R3 
R3*  < OM ( I - 1 ) +3 . *0M  < I ) )*R3 
G 1 ( I ) =P1 +G+R1 
G2( I )=P2+R2 
G3< I )=P3-Q+R3 

CU<  X )«-Pl*U(I  +  l >-P2*U( I )-P3*U( 1-1 ) 

C  THE  DIFFUSION  TERM 

AU< I )=2./0MD 

RU<  n=SC<  1-1  )*AU<  I  )/<0M<  I  )  -OM  <  I  -  1  )  ) 

AU< I )=SC< I >*AU( I )/<0M( 1  +  1 )-0M< I )  ) 

IF  <  NEQ . EQ . 1 )  GO  TO  20 
DO  15  J= 1 . NPH 

C(.l.  I  )'=-Pl*F<  J.I  +  1  )-P2*F(  J.I  )-P3»F(  J.I-1  ) 

CALL  SOURCE ( J » I . CS  »  D ( J  » I ) > 

C<  J. I )=-C<  J. I )+CS-F<  J. I )*D< J. I > 

A< J, I )=AU( I )/PREF< J) 

B( J. I )=BU< I >/PREF< J> 

15  CONTINUE 


C  SOURCE  TERM  FOR  VELOCITY  EQUATION 

:  o  si  ( i  >*ofdx*dx 

s:1  <  I  '  -P24S1  (  I  >/<RHO<  I  )  *U  (  I  )  ) 
S3<I>=R3*S1<I ) / ( RHO  < I-i >*U< 1-1 ) ) 

SI < 1 >«P1*S1 < I  )  /  (  RHO  (  IM  >*U<  141  >  ) 

C  *J  <  I  )*-ClM  I  >  -  2  .  *  (  S 1  (  I  >  ♦  S  2  <  I  >  ♦  S  3  <  I  >  ' 
SI < I >*S1 ( I )/U( I ♦  1  » 

S2< I >  *S2  < I )/U( I  ) 

S3 ( I >=S3< I >/U<  1-1  ) 

25  CONTINUE 

C  COEFFICIENTS  IN  THE  FINAL  FORM 

DO  30  1*3. NP 1 

Rl  *■  1  .  /  <G2<  I  )  4AU  <  I  >  -f  fiU <  I  >  - S 2  <  I  )  ) 

AU< I )  =  < AU< I >4S1 ( I )-Gl < I >  >tRl 
BU< I >*<BU< I >4S3< I ) -G3  < I > )*RL 
30  CU( I >=CU< I ) *RL 

IF(NEQ.EQ.l)  GO  TO  40 
DO  35  J*1,NPH 
DO  35  I *3  *  NR  1 

RL*1 ./<G2< I  >4A< Jt I )♦»( J» I )-D< J. I ) > 
A(J.I)*<A< J. I ) -Gl < I ) >*RL 
B<J.I)*<B<J.I)-G3<I>>*RL 
C<  J. I )*C<  J. I )*RL 
35  CONTINUE 
40  CALL  SLIF 
RETURN 
END 


SUBROUTINE  VEFF < I » I P 1 . EMU ) 

C0MM0N/BLK2/N » NP1 . NP2.NP3.NEQ.NPH.KEX. KIN.KASE » KRAD » PE  I » AM  I » AME  » 

1  DPDX . XU. XD » XP. XL .OX  » INTG.CSALFA .HGS. CINF  < 10)  »PREF< 10) . 

2  PR( 10) ,P< 10) »U<050) *F< 10.050)  »R(050) » RHO (050) . OM < 050 > . Y < 050 ) . 

3  DEN(OSO) »AMU<050) 

C0MM0N/BLK3/KRP » MACH » KR » INDI ( 10) . I  NOE  < 10) . SUMR » AK . ALMG . BET A . TAU I 

1  .TAUE.VINF.DELTA.TU.TINF.QS. YL. UMAX. UMIN.FR. YIP. YEM.TR.HW. 

2  HF< 10) »  GAMA ( 1 0 ) » A Jl < 1 0 ) » A JE ( 1 0  >  » SC (050 ) . AU(050) »BU<050) .CU(050)  » 

3  A<  10.050)  »B<  10.050)  »C<  10.050)  ,TEMF'E<  050)  .  TEMP  (  050  )» PO  (  050  )  . 

4  AMACH(OSO) .CP (050) 

COMMON /BLK4/NH » NNO »  NNP . NO  »  NN .N02.NN2.NE. IND< 1 0 ) . STO . AKS . RT . FT . 

1  AMT . DRAGC . STNO . TAU I U . QU . UGU . UGD 
C  THIS  SUBROUTINE  USES  THE  MIXING-LENGT  HYPOTHESIS 

Al.  =  ALMG*YL 
YM=  (  Y  (  1  )+Y<  IF-1  )  )*.5 
TFTYM.LT  ,AI./AK)AL  =  AK*YM 

10  EMU*  .  5*  <  RHO  (  I )  4RHG <  IF' 1 )  )  *AL*AL*ABS(  (  U (  I ) -U (  IP  1  )  )  /  (  Y  (  I  )  -  Y  <  I P 1  )  )  ) 
RETURN 
END 
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■  mm  mi  t  ?  mi  snni.  i  f  i  .1 .  i .  t  s .  ps  > 

I  'It.  rONM  fVAUON  (If  f  T  AON  A  T  I  ON  FNTHAIFt 

r.ViTKlN  list  C  f)  N  ‘‘  1  ^  T  F  NT  UNITS 

!M(  P'M  »  f  nim(  T  or  I  WITH  t  IS  NFGLFCTIP 

|  ihhun  hi  M  -  I  r-t  .  IS.AS.FS.  1 £  MPA  (  050  )  .  RMOA  <  050  >  •  T  A  <  050  »  .  i  ►.’  <  050  '  . 

1  f  !•:*  >  Or.n  t  .  iau:m  050  >  .  T  AIJ3<  050  )  .  7fc?  •  050  >  .  ?N  3<  050  > 

L  |l«NHN,  |<|  K  ,'/N.  Mf  1  .Nl  I’.trr  J.NIO.MI  H.KfK.MN.KASI  .  M>  Ap.  »  f  J  .  AH  1  .  AHf  , 
1  [If  (■»  .  XII.  nil.  yf  .  VI  *  fly  ,  INTO  .  CSAI  F  A  .  MGS  *  C  I  NF  <  I  0  >  *  PRFf  <  1  0  t  . 

'  PR  (  1  0  )  .f  (  I  0  >  .  U<  050  »  .f  <  10.050  >  .  R  <  050  )  .RMO<  050  >  •  QH<  050  )  .  V  <  050  >  . 

3  t>FN<  05«  >  •  AMU  <  050  * 

COMMON/ K  K3  /*  Rf‘  .  MACM  »hhi  I  MM  <  10)  .  INpE  (  10)  . SUNK. Ak.Al  MG  .  PF  T  A  .  T  AH  I 

1  .  T  AUF  .  V  INF  •  I*  F  I  T  A  .  T  W  .  T  INF  #  OS  *  Y  L  .  UMAK  ,  UNINiFhi  V  I  F'  »  Y  f  M  »  TP , MU . 

2  MF  <  1  0  )  .  GAMA  <  10)  •  A  .1 1  «  lOl.A.tf  (  1  0  )  *  SC  <  050  )  .  AU  <  050  >  » PU  <  050  )  •  C(J «  050  )  . 
J  A  » 1 0 .050  > . p<  I  0. 050 ) .( < 1 0.050  > . TCHPF  < 050 > . TFMF  <  050  > . f  0 <  050 ) . 

A  AMA(  H  <  050  '  .  Cf- «  050  > 

roMMON/M  k 4 /NM , NNO . NNf  •  NO  *  NN  .  NO!’  ♦  NN2  *  ME  *  INI"  10)  «STO.AKS»RT  .FT  . 

1  AMT  .  PPAGC .STNO.  T  AU  I  U  .  OW  »  UGU  .  UC.P 
M  Ml’ NS  I  ON  kCS<  200) 

RF  N-0 . 0  7  5 

I f  <  J . MF . 1 >  r 0  TO  20 

1'  S  r  SC  <  I  >*(()<  1  ♦  1  )  *  U  <  I  4  1  )-U<  I  >  *  U  (  I  )  )  /  (  OH  <  I  ♦  1  )  -  OM  (  I  )  ) 

CS-CS-SC  < I  -  1  > • <U<  I  >  *U<  I > -U<  I  -  1  >*U< I  -  1 >  >/<  0H< I > -0H(  I  -  1  >  > 

CS»  (  1  .  -  l  .  /PREF  (J|)*CS/<0H<MI>-0H<I-1)> 

IF  < KF . FO. T )  GO  TO  15 
A  S  "  0 . 0 
PS  ~  0 . 0 

if  <  NR  ,  T  (),  1  )  CAL  L  RAPCOM  (  I  ) 

I'S*  AS4CR 

I'Sr  pc 

(.0  TO  25 
CUNT  1 NUt 

7h  =  11 5 . * • RHO< I  >/PEN)*<  TEMPT I )/( 10. **4.  ) >**5. 

RCS< I >  =  2. » 1  . 71 4E-9»ZK*<0S4TEMP< 1 >**4. 4TEMP(NP3  >  **4 . -2*<  TEMP< I  > 

1  *44  .  )  ) 

CS-  .  21  62*32. 24 PCS < I ) / < RHO< I >*U< I ) >  4CS 

PS*  <  2. *  1 . 714E-9*<  0. 0575*<  TEMPT  I  >/ < 10. **4 .  >  >**4 . *  <  TEMF  < 1 >**4 . 4 

1  TEMP < NP3  >*»4 . -2. *<  TEMPT  I  )**4  .  )  ) / PEN-B . * <  ZK /RHO <  I>>*TEMP<I)**3.)/ 

2  <CP< I )*U( I ) ) >*.2162*32.2 
IF(  I  .E0.NP1 )  GO  TO  30 

GO  TO  25 
CONTINUE 
CS^O.O 

ns=o.o 

CONTINUE 
RETURN 
CONTI NUf 

/M-l  1  5  .  *  <  R  H  0  (  1  )/PEN>*< TEMPT  1)/(10.**4.))**5. 

ZKNF'3  =  1  15.  *  <  RHOT  NF'3)/PEN)*  <  TEMP  T  NP3 ) / T  10.**4.))**5. 

PC  1 *2 . *  1 .  7  14E-9*ZK1 <<TEMP<NP3 >**4. -TEMPT 1 )**4. ) 

RCNf'3  =  2  .  <  1 .714E-9*ZKNP3*<  TEMP  <  1  )  **4  .  -TEMP  <  NP3  )  **4  .  ) 

SUMR*0. 0 
DO  35  J-4.NP1 

SUMR  *  <  RCS  < I ) +RCS  T I-1))*<Y<I)-Y(I-1))*. 5  +  SUMR 
CONTINUF 

SUMR*SUMR4(RCST3)4RC1 )*<Y<3)-Y<1> ) * . 54 < RCNP3 4RCS < NP 1 ) )*<Y(NP3>- 
1  Y<NP1))*.5 

SUMR»TGS4.5*SUMR41 .714E-9*TEMF(1)**4. )/3600. 

RETURN 

ENI.  82 


inn."u  r  r  m(  st  r  f 

IMMT  N/ p|  K2/N. NF  1  .NF  ?.NF  3.NE Q.NF  H.Kf  X . K  IN.KASC  .KF Ap.F f  T . AMI  . AM  I  . 
i  pr  in  .  xu.  xp,  xf  .  xt  . px . intg. c sai r a.mgs . c inf  <  i o  >  .  pfef  <  i  o  > . 

ff< io > . f  < i o > . u <  rtf.o  > .  r  <  i  o .  nso  > .  f  <  050  >  .khui  oso  > .  om  <  oso  >  •  v  <  oso  >  . 

3  Pf  N<  050  )  .  AMU<  050  ) 

COMMON /Ht  K  3/KPF-.MACM.KF  .  INPI  (  10)  .  1  NPE  (  10)  .  SUM*  .AK.AlMG.PfTA.TAUI 
1  .  TAtlf  .  VlNF  .  PF1  TA.TW.TINF.OS.VL.  UMAX  . UM I N . r F . Y 1 F . YE M .  TF  .  HU  . 

?  MF  »  10  >  .GAMAi  10  >  .  AJ1  <  1  0  >  .  A.if  <  1  0  )  .  SC  <  050  >  .  AU<  050  >  .  PU<  050  >  . f  U<  050  >  . 

3  AM  0.050  )  .  P  (  1  0 . 050  )  .C  <  10. 050  '  .  Tf  MFC  <  050  >  .  Tf  Mf(  050  )  .  F  0<  050  )  . 

4  AMACHI 050 ) .CF < 050 ) 

C  0MM0N/PLK4 /NH . NNO . NNF  . NO . NN , NO? . NN2 . Nf , I NP <  1 0  >  *  S T 0  *  AK S . R T , F  T , 

I  AMT  ,  pFAGf  .  ST  NO  ,  T  AU  I  W  ,  OU  .  UGU  .  UG  P 

r  St  IF-  COEFFICIENTS  NEAR  THE  I  POUNPARY  F  OF  VELOCITY  EQUATION  EQUATION 
riii  ?  >  *o . 

CU«  NP?  >  *  0 . 

GO  TO  ( 10. 15.20)  .KIN 
10  PU ( 2  > -  0 . 

All  (  2  )  **  1  .  /  <  1  .  ♦  2  .  •  BE  T  A  ) 

GO  TO  30 

15  S0'-84*U(1)*U<1)-12.»U<1)*U(3>+9.*U(3)*U(3) 

PU(2)*8.*<2.tU< 1 ) ♦ U  <  3  > > / <  2 . *U  < 1 ) 4 7 . *U ( 3 ) 4 SORT < SO > > 

AU<  2)  =  1 . -PU(2> 

GO  TO  30 
20  PU <  2  >  r0 . 

CALL  VEFF (2.3. EMU  ) 

AK  1  - 1  .  /PX  -PPDX/  <  RHO  <  1  >*U< I >*U< 1  >  > 

AK2  =  -U(  1  >*AK14PPDX/<RH0<  1)*U(  1  )  ) 

AJ*F?HO<  1  )tU(  l  >*.25*<Y(2)4Y<3>  >**2/EMU 
IF(KFAP.EQ.O)  GO  TO  25 
AU<2)*2./<2.4AJ*AK1  ) 

CU<2)=-.5*AJ*AK2*AU<2) 

GO  TO  30 

25  CU ( 2 ) s 1 ./(2.43.BAJ4AK1) 

AU(2)*CU<2)*<2.-AJ*AK1 > 

CU(2)*-CU(2)*4.»AJ*AK2 

C  SLIP  COEFFICIENTS  NEAR  THE  E  BOUNDARY  FOR  VELOCITY  EQUATION 
30  GO  TO  (35.40.45) »KEX 
35  AU  (  NF‘2  )  *0  . 

PU(NP2)*1./<1.+2.*PETA) 

GO  TO  50 

40  SQ*84.*U<NP3)*U(NP3)>12.*U(NP3)*U(NP1)+9.*U<NP1 >*U<NP1  ) 

AU(NP2)=8.*(2.*U(NP3)+U(NP1 ) ) / ( 2 . *U ( NP3 ) 4 7 . *U ( NP 1 )4SQRT<SQ> ) 

PU ( NP2  >  =  1.-AU(NP2) 

GO  TO  50 
45  A  Li  (  NP2  )  =0  . 

CALL  VEFF  <  NP 1 »  NP2  »  EMU ) 

PK1=1./DX-DPPX/<RH0<NP3)*U<NP3)*U(NP3> ) 
PK2=-U(NP3)*PK1+DPDX/(RH0(NP3)*U(NP3) ) 

BJ=RH0(NP3)*U<NP3)*.25*<2. *Y < NP3 > - Y < NP 1 )-Y(NP2) )** 2/EMU 
CU(NF'2)~1  ./(2.4-3.PBJ4BK1  ) 

BU<NP2)=CU<NP2)*<2.-BJ*BK1 > 

CU(NP2)=-CU(NP2)*4 . *BJ*BK2 
50  IF(NEQ.EQ.l ) RETURN 
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C  SLIP  COEFFICIENTS  NEAR  THE  I  BOUNPAR Y  FOR  OTHER  EQUATIONS 

HM  1 0*  J-ltNFH 
i.i  J,  ?  )  0  . 
ft  l,NP2>=0. 

SO  10  (55. AS. 70) ,RIN 
CALI  FBC  <  xit,  J.  INDI  <  J)  »QI> 

I4MNPI(J>.EQ.1>  GO  TO  60 
A  J 1 ( J)*OI 
A<  .1.  2)  =1  . 

B( J,2>=0. 

n<J.?>*8. *<1.42.* BETA)  *PREF< J)*AJI  <  J > / ( AK*AK*BETA* < 1 . 4BETA  > * ( 1 . 4 
1  BE T A ) * ( 3 . *RHO < 2>4RH0<3> )*U<3)  ) 

GO  TO  80 
AO  F ( J  » 1 >  *GI 

A( J»2)*<  1  . 4  BE  T  A-GAHA  <  J ) )/< 1 . 4  BE  T  A4GAMA ( J )  ) 

B< J.2)*l .-A( J.2) 

GO  TO  80 

65  A<J.2)*(U<2)4U<3)-8.*U( 1 ) )/<5.*(U(2)4U(3) )48.*U<1 ) ) 

GF»< 1 ,-F  REF< J) ) / < 1 . 4PREF  <  J ) ) 

A< J»2)«<A<  J.2)4GF)/< 1 .4A< J.2)*GF ) 

B< J.2>«1 .-A< J.2) 

GO  TO  80 
70  B ( J »  2 ) =0 . 

CALL  SOURCE (J.l.CS.BS) 

AK1*1 ./DX-DS 
AK2*-AK1*F< J. 1 )-CS 
A JF*A J*PREF  <  J ) 

IF(KRAD.EQ.O)  GO  TO  75 
A< J.2)*2./<2.4AJF*AK1 > 

C< J.2>=-.5*AJF*AK2*A< J.2) 

GO  TO  80 

75  C<J»2)=1./<2,43.*AJF*AK1) 

A< J,2)=C< J.2)*<2.-AJF*AK1 ) 

C< J.2)*-C< J»2)*4.*AJF*AK2 

C  SLIP  COEFFICIENTS  NEAR  THE  E  BOUNDARY  FOR  OTHER  EQUATIONS 

80  GO  TO  (85.95. 100) »KEX 
85  CALL  FBC(XD.J. INDE< J) . QE) 

I F ( I NDE <  J ) . EQ • 1 )  GO  TO  90 
A JE ( J ) =QE 
B< J.NP2)*1  . 

A ( J , NP2 ) *0 . 

C<  J»NP2)=-8.  *(  1.42.*BETA)*F'REF(J)*AJE(J)/(AK*AK*BETA*(1 . 4  BET  A  )  * 

1  ( 1 . 4BETA)*(RH0(NP1 ) 43 . *RHO < NP2 ) )*U<NP1 ) ) 

GO  TO  105 
90  F<J,NP3)=QF 

B<  J.NP2)=< 1 .4BETA-GAHA< J) ) / < 1 . 4BETA4GAMA < J )  ) 

A< J.NP2)=1 .-B< J.NP2) 

GO  TO  105 

95  B  <  J » NP2  )  =  (  U  <  NP2  )4U(NP1  )-8.#U(NF'3)  )  /  (  5  .  *  <  U  ( NP2  )  4U  (  NP 1  )  >48.* 

1  U  <  NP3 ) ) 

GF  =  <  1 .  -F'REF  <  J  )  )  /  (  1 . 4F‘REF  (  J )  ) 

B  <  J  »  NF'2 )  =  <  B  <  J ,  NP2  )  4GF  )/<1.4B(J»  NP2  )  *GF  ) 

A  <  J » NP2 )  =  1 . -B <  J . NP2 ) 
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GO  TO  105  . 

100  A ( J  f  NP2 ) =0 • 

CALL  SOURCE< J.NP3.CS.DS) 

BK1  =  1 ,/DX-DS 
BK2=-BK1*F(J»NP3)-CS 
BJF=BJ*PREF< J) 

C(  J,NP2)  =  1  . /<2.+3.*BJF*BKl  > 

B<  J»NP2)=C<  J»NP2)*<2.  -BJF#BN1  ) 
C(  J»NP2)=-C<  J»NP2)*4.*BJF*BK2 

105  CONTINUE 
RETURN 
END 


SUBROUTINE  FBC ( X , J , IND r A JFS ) 

COMMON /BLK 3 /KRP »  MACH » KR » INDI ( 10 ) » I NDE < 1 0 > . SUMR . AK . ALMG » BETA . T AUI 

1  .TAUE.VINF.DELTA.TU.TINF.OS.YL.UMAX.UMIN.FR.YIP.YEM.TR.hu. 

2  HF < 10) .GAMA < 10) .AJI < 10) . AJE( 10) .SC (050) . AU ( 050 ) . BU ( 050 ) .CU(050> , 

3  A ( 10.050) »B( 10.050) . C ( 10 . 050 ) » TEMPE ( 050 ) . TEMP ( 050 ) .P0(050) . 

4  AMACH(050) »CP<050) 

IND  =  1 

C  H  MUST  HAVE  UNITS  FT . FT/SEC . SEC 

IF  (J.EO.l)  GO  TO  105 
GO  TO  110 
105  A JFS  =  HU 
110  IF ( J ♦ EQ . 2 )  GO  TO  115 
GO  TO  120 
115  A JFS= 1 .  O-CU 
120  IF ( J . EQ . 3 ) GO  TO  125 
GO  TO  130 
125  A JFS=CU 
130  RETURN 
END 
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non 


SUBROUTINE  EXPNI < RES . N f X , I ER > 

THIS  SUBROUTINE  EVALUATES  THE  EXPONENTIAL  INTEGRAL  FUNCTIONS  IN 
THE  RADIATION  SOURCE  TERHS 
TEST  OF  RANGE 
I  ER=0 

IF(ABSCX) .LE. .00001 )  GO  TO  10 
GO  TO  25 

10  IF ( X ) 1 5  »  20  f  20 
15  X=-. 00001 

GO  TO  25 
20  X=. 00001 

25  IFLX-4.)  30 »  55  *  35 
30  IF (  X  +  4 . )  65»55»50 

C  ARGUMENT  IS  GREATER  THAN  4 

35  ARGM./X 

RES=EXP ( -X ) * ( ( <((((( . 000944276 144ARG-. 0049362007 )*ARG+ 

1  .011 723273 )*ARG-. 01 7555779 >*ARG+. 02041 2099 )*ARG-. 022951 979 )tARG+ 

2  .031208561 J4ARG-. 062498588 )*ARG+. 24999999 >*ARG 
40  IF(N.EO.l)  RETURN 

NM1 =N-1 
EXPX=EXP  < -X ) 

XN 1  *  1  . 

HO  45  N1  =  1 f  NM 1 
RES=(EXPX-X*RES)/XN1 
45  XN1 ~XN1  +  1  . 

RETURN 

C  ARGUMENT  IS  ABSOLUTELY  LESS  OR  EQUAL  4 

50  IF ( X ) 55  f  60  »  55 

55  RES=-ALOG(ABS(X) )-( ( (<<<<<<<((( . 1 03 1 7602E- 1 1 *X- . 157986/ 5E- 1 0 > *X+ 

1  . 16826592E-9)*X-.21915699E-8)*X+. 27635830E-7 > *X- . 30726221E-6 > *X+ 

2  .30996040E-5)*X-. 28337590E-4 > *X+ . 231 48392E-3 ) *X- . 001 6666906 >*X  + 

3  .01 0416662>*X-. 055555520) *X+. 25) *X-1 .0)*X-.57721 566 
GO  TO  40 

60  RES= 1 . E37 

GO  TO  40 

C  ARGUMENT  IS  LESS  THAN  -4 
65  I E  R  = 1 
RETURN 
END 
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SUBROUTINE  SOLVE  <  A  *  B  .  C  *  F  » NF'  3  ) 

THIS  SOLVES  EQUATIONS  OF  THE  FORM 
F  < I )^A( I >*F( 141 >  +  B( I ) *F ( I  -  1 ) EC( I  ) 

FOR  I=-?,NP? 

HI  HENS  I  (IN  A  (  NF'3  )  f  P(  NF'3  >  ♦  C  <  NF'3  )  »  F  <  NF'3  ) 
NF2=NP3-1 

B<  ?>  -  B  (  2 ) *F (  1  )  +  C  (  2  ) 

0(1  10  I  =  3  t  NR  2 
T  =  1  .  /(  1  .  -  B  (  I  )  *  A  (  I  -  1  )  ) 

A (  I  )  =  A(  I  )*T 

B( I (B( I )*B< 1-1  )  +  C ( I ) ) *T 
DO  15  1=2. NR? 

J=NR2- 1+2 

F(  J)=A(  J)*F(  J-fl >FB( J) 

RETURN 

ENH 


REAL  FUNCTION  EBLT(K.UV) 

REAL  UV ( K ) 

THIS  FUNCTION  SUBPROGRAM  CALCULATES  PLANCKS  RAO I AT  I  ON  FUNCTION 
COMMON/ BLK 1/LST.TS. AS  *  BS . TEMP A ( 050 ) »RH0A(050) »YA(050)  »EB2(050)  . 

1  EB3(050) »TAU2(050) »T  AU3 ( 050 ) »  ZK2 ( 050 ) »ZK3(050) 

COMMON /BLK3/KRF'  ,MACH»KR.INDI( 10) , I NBE ( 10) .SUMR.AK.ALMG.BETA.TAUI 

1  .TAUE.VINF.OELTA.TW.TINF.OS. YL.UMAX.UMIN.FR, YIP, YEM.TR.HUt 

2  HF< 10) .GAMA< 10) , AJI (10) .AJE( 10) .SC (050  > , AU(050 ) . BU( 050) , 

3  CU(050) *A( 10. 050) ,B(10,050) »C( 10,50) .TEMFE (050) .TEMP (050) . 

4  POT  050  )  » AMACH  (  050  )  » CF‘  (  050  ) 

C 1 = 1 . 1 87E  08 

C2=25896 « 0 

IF  (LST.EQ.O)  GO  TO  10 
T=TEMPA(LST) 

GO  TO  15 
T  =  TP 

CONTINUE 

ARGA=C2/(UV( 1 )*T ) 

IE  (ARGA.GT.50. )  GO  TO  20 

BAH  =  Cl/( (UV( 1 )t*5. ) ♦ ( EXP ( ARG A ) - 1  .  ) ) 

EBLT  =  BAD 

RETURN 

CONTINUE 

BA1iEX  =  AL0G(C1  )-(5.tAL0G(WV(  1  )  )  +  ARGA  ) 

BAB=EXP(BABEX) 

EBLT=BAH 

RETURN 

END 


»  t'NC  »  I  ON  *.HUnU  (  I  ,  .1  > 
kf  rtl  *1  ' 

I  ennOW  'APf  A  1  'H.  nil  <  8  A  <  1  O  > 
MHfWSION  Al<?0' 

£  *  I  F  PNAl  PN 
HIM  1  )=0.06943184£0 
net  7)-0 , 33000984E0 
ne ( 3  >  =0 . 66999052F0 
h" ( 4 >  =0. 93056816F0 
AIM)*?. 0053 
A1  (2)  =  1 .53831 
A  1  (  3  )  *0 . 6946*; 

A  1  <4>=0. 23156 
A  1  <5)«0. 04996 
A  1 (6)*0. 00713 
A 1  <  7 ) *0 . 00089 
A1 ( 8  )=0. 0 
I'O  10  I  P=  1  »  4 
nu<  IP-M  )  *-HU<  IP  ) 

N^P 
I  I  I  =  I 

IF  (  I  .LT.O)  III*-Ii4 
JJJ=  J 

IF ( J.LT.O) JJJ*-J44 
S*1 .0 

l<0  20  H*  1  *  N 

nnn  =  n 

S*S+A(M)*PN(III» JJJ) 

snunu=s 

RETURN 

FND 


FUNCTION  PN(ItJ) 

COMMON/A RFA1/M»MU (8) *  A  ( 10) 

EXTERNAL  PNCOS 

PXRC=DMLIN<PNC0S»0.0»2*3. 1 4159* 1 » 300 *  0 . 0  *  0. 001 » IER > 

PN* 1 . / ( 2 . 0*3 . 1 4159 ) *PXRC 

RETURN 

ENH 


FUNCTION  FNCOS<Kf PHI  ) 

C OHM ON /ARE  A 1/M*MU(8)  . A< 10) 

REAL  HU  »  PHI ( K  ) 

FACT<Z >=EXP< -2 ) *Z**Z*SQRT<  2. *3 . 14159*Z  > 

COSTH( APHI .  X  » Y )=X*Y+< 1 . -X**2 ) **0 . 50* < 1 . -Y**2)**0.50* 

J  COS(APHI) 

PNC=0. 0 
N- M 

N  P  1  =  N  +  1 
ZN=FLOAT(N> 

I«0  15  MM  =  1 » NP1 
ZM=(HM-1 )/2. 

PNC= ( - 1 . )**ZM*FACT (ZN)/(FACT (ZM)*FACT (ZN-ZH) ) *FACT ( 2 . *ZN-2 . *ZM ) 
1  <FACT<ZN)*FACT<ZN-2.*ZM) ) *COSTH < PHI < 1 ) . MU < I ) . MU ( J ) ) ** ( ZN- 2 . * ZM ) 
PNC0S= ( 1 . / ( 2 . **M ) ) *PNC 
RETURN 
END 


SUBROUTINE  RADF'RM 

COMMON/BLKl/LST •TStASrBS»TEMPA(050) *RH0A(050) r YA(050 >  ,EB2<  050  ’  . 
1  EB3 ( 050 ) r  TAU2  <  050 ) . TAU3 < 050 ) r ZK2 ( 050 > » ZK3(050) 

C0MM0N/BLK2/N.NP1  .NP2f NP3»NEGf NPHtKEX. KINt KASErKRAD.FE I  ,  AMI  .  A«I 

1  DPDX>XU>XD»XPt XL»DX* INTG . CSALFA f HGS » C I NF < 10) tPREF (10), 

2  PR< 10) »P< 10) ,U<050) ,F< 10 » 050) ,R(050) .RH0(050> *0M< 050 > . >'  0?  2  • 

3  DEN(050) »AMU(050) 

C0MM0N/BLK3/KRP»MACH»KRf INDI ( 10 ) » I NDE ( 1 0 ) t SUMR , AR . AL  HO • BC  ’ a  •  ’ 

1  . TAUE. V INF .DELTA, TU. TINF.GSr YL  . UMAX  rUMIN.FR, YIP.YEH.TF.mw. 

2  HF ( 10) *  GAMA ( 10) . AJI < 10) » AJE( 1 0 ) .  SC ( 050  > »  AU ( 050  )  . BU <  ?r •  <.'  • 

3  CU<050) »A<10»050) » 10»050) »C( 10? 50) * TEHPE < 050 ) . TF 

4  P0(050) »AMACH(050) . CP  <  050 ) 

COMMON/ BLK5/KR I» P INF  r  B IFF  <  50  > f NOPT.NSKAT 
COMMON/SAN J/SIGMEXfSIGMABfPDEN, SIZE 
DIMENSION  AUX<  300) r TRAN2  <  50  > r  TRAN3<  50  > 

DIMENSION  EXTN2<50) »EXTN3(50) »AK2<50) . AK3( 50 ' . 7RAT > *  •  ’  *  *  •  ' 

1  EX2 ( 50 ) » EX3 (50 ) » AR2 ( 50 ) f  AR3  <  50 ) f  AD2 ( 50 ) . AD  3  < 1 ;  . 

2  BD2(50) f  BD3<50 ) f  ZT2( 50 ) . ZT3<  50 ) »  AA(  4  > . YO  -  *  o 
EXTERNAL  EBLT 

EXTERNAL  SMUMU 
EXTERNAL  AKP2 
EXTERNAL  AKP3 
BEN=0.075 

I F  (  NSK AT  .  F.Q  .  1  >  GO  TO  3 
I F ( NSKAT . EQ . 3  >  GO  TO  2 
AA ( 1 ) = . 1 739274 
AA ( 2 ) = . 326072* 

A  A  (  3  )  -  .  3  2  6  0  7  2  A 
A  A  (  4  )  =■  ,  j  739274 
DO  1  1-1.4 

SCN-0 , 0 
DO  1  J -  1 . 4 
SEN- SCN4AAM  >*<  AA< 

GO  TO  3 


MICROCOPY  RESOLUTION  TEST  CHART 


NJ  r  J 


?  SCN-0.5 

3  DO  10  L  =  3  *  NP1 

TEMPA<  L-l )=TEMP  <  L ) 

RHOA (L-l )=RHO(L) 

YA(L-1 )=Y<L> 

YD<L)=Y<L+1 )-Y<L) 

10  CONTINUE 

YD(2)=(Y<3)-Y(1> >/2.0 
YD ( 1 )  =  Y D  <  2 ) 

TEMPA< 1 ) = TEMPI  1 ) 

RHOA  < 1 ) =RHO  < 1 ) 

YA ( 1 ) =Y ( 1 ) 

TEMPA(NP1 )=TEMP<NP3> 

RHOA  <  NP1 ) =RHQ ( NP3 ) 

YA  <  NP1 ) =Y<  NP3 ) 

T  AU2  < 1 ) =0 . 0 
T  AU3 ( 1 ) *0 * 0 
DO  25  L-l t NP1 

ZK2(L >=<4. 37/0.03281 )*< (RHOA(L) /BEN) **1 .009)* < (TEMPA(L)* 

1  0,00018)**2.85) 

ZK3(L)  =  ( .04985/0.03281 )*<  < RHOA < L > /BEN )** i .205)* ( (TEHPA(L)* 
1  0 . 00018 ) **5 . 47 ) 

LST=L 

VFRA=F(3,L)*RH0(L)/PDEN 
NPART=VFRA*YD(L ) /SIZE 
EXTN=SIGMEX*NPART+SIGMAB*VFRA 
SKA  =SIGHEX*NPART 

EB2  <L ) =DMLIN  <  EBLT  t .  05  » .135»1»300»0.0»0.001»IER) 
EB3(L)=DMLIN(EBLT» . 135 » 20 . » 1 » 300 r 0 . 0 » 0 . 001 r IER ) 

ZKA2  <  L )  =  SIGMAB*VFRA  +  ZK2  <  L ) *  < 1 ♦ O-VFRA ) 

ZKA3  <  L )  =  SIGMAB*VFRA+ZK3(L)*<1. O-VFRA) 

IF < NSKAT . EQ . 1 )  GO  TO  12 
EXTN2(L)=EXTN+ZK2<L>*< l. O-VFRA) 

EXTN3 < L )  *  EXTN+ZK3<L>*<1. O-VFRA) 
AK2<L)=<ZKA2<L)*(ZKA2<L)+2.0*SCN*SKA) )**0.5 
AK3(L)=(ZKA3(L)*(ZKA3(L)+2.0*SCN*SKA))**0.5 
GO  TO  15 

12  AK2(L)=ZKA2<L> 

AK3 ( L)=ZKA3<L) 

15  IF(L.EQ.l)  GO  TO  25 

TAU2  <  L )  =  <ZK2(L)+ZK2<L-1) )*< YA(L)-YA(L-1 ) ) *0 . 5+TAU2 < L- 1 > 
TAU3<L  >  =  (ZK3(L)+ZK3(L-1 ) )*< YA(L)-YA(L-1 ) ) *0 . 5  +  TAU3 < L- 1 ) 

0  CONTINUE 
5  CONTINUE 

EX2<1)=(EB2<1)+EB2<2> )/2.0 
EX3< 1)*<EB3<1)+EB3<2) )/2.0 
AR2(1)=EB2(1)*2.0*ZKA2(1) 

DO  30  1.-1  »NP1 
AR2(L)«EB2(L)*2.0*ZKA2<L) 

AR3<L)*2,0*ZKA3<L)*EB3<L> 

30  CONTINUE 


90 


00  70  L  =  1  » NP1 

AD2<L)=<2. 4AKP2  < - 1 .  » 0 . » - Y ( L  + 1 > »AK2(L) »ZKA2<L> > #EX2 < L > -ZKA2 < L ) * 

1  AR2(L)/<  AK2  <L) )**2*( AKP2 ( -1 • »0. t-Y<L  +  l ) , AK2 ( L ) , ZKA2 < L ) > -AKP2 < 1 . » 

2  0.  t-Y(L)  »AK2(L)  »ZKA2(L)  )  )  )/(AKF‘2(-l  .  »  1  .  ,  Y  (  L  ) -Y  <  L  + 1 )  »  AK2  ( L  )  » 

3  ZKA2<L>  >-AKP2< 1 . t\ . » Y < L  +  l > - Y < L >  * AK2 < L >  » ZKA2 < L >  > ) 

AD3 ( L )  =  <  2 . *AKP3 ( -1 . ,0. » - Y < L  + 1 ) » AK3 ( L ) * ZK A3 < L ) ) *EX3 ( L ) -ZKA3 ( L ) * 

1  AR3  ( L ) / <  AK3 ( L ) ) **2* ( AKP3 < -1 . » 0 . , - Y ( L  + 1 ) r AK3 < L ) » ZKA3 < L ) ) -AKP3 ( 1 . , 

2  0. f-Y<L>  f AK3(L) #ZKA3(L) > ) )/<AKP3(-l . , 1 . , Y < L ) -Y ( L+l ) »AK3(L) » 

3  ZKA3(L> )-AKP3< 1 • » 1 « ♦  Y  <L  +  1)-Y(L) »AK3(L) »ZKA3<L) ) ) 

BD2  <L)=(2«*EX2(L) -AD2  <L)*AKP2(-1 . »0. » Y ( L ) » AK2 ( L ) » ZKA2 ( L ) ) - 
1  ZKA2 (L>#AR2(L)/<AK2<L>  >**2»0) /AKP2 ( 1 ,t0, t -Y < L ) » AK2 ( L ) » ZKA2 ( L ) ) 
BD3(L)  =  (2.*EX3(L)-AD3<L)*AKP3<-1 .>0, t Y ( L ) » AK3 < L >  » ZKA3 ( L ) ) - 
1  ZKA3<L)*AR3<L)/<  AK3  <  L ) ) **2 . 0 ) /AKP3 < 1 . »0. »-Y(L) »AK3(L)»ZKA3(L) ) 

32  IF(EX2(L) . EQ . 0 . 0 )  GO  TO  35 

ZT2<L)  =  1 .0-<0.5#AD2<L)*<AKP2<-l . *0. »Y<L+1 ) *  AK2 ( L ) »  ZKA2 ( L ) ) - 

1  AKP2 ( 1 . » 0 ♦ »Y(L)fAK2(L) »ZKA2(L>  >  >+BD2<L)*0.5*<  AKP2  <l«»0»r-Y(L  +  l) 

2  rAK2<L)»ZKA2(L))-AKP2<-l . » 0 . » -Y < L ) » AK2 < L ) > ZKA2 < L ) ) ) >/EX2(L) 
TRAN2<L)  =  < AD2(L)*AKP2<~1 »»0«»Y<L41)»AK2<L>  »ZKA2<L) >  + 

1  002 < L ) *AKP2 ( 1 . ? 0 . »-Y<L  +  l )t AK2(L) tZKA2(L)  )  + 

2  ZKA2(L)*AR2<L)/<AK2(L)**2> )/<2»*EX2<L> ) 

GO  TO  40 

35  ZT2  <  L ) =0 . 0 

TRAN2<L  )  =  1 .0 

40  IF<EX3(L)  ♦EG.O.OGO  TO  45 

ZT3(L)=1.0-(0.5*AD3<L)*(AKP3(-1 . r 0 . » Y ( L+l ) f AK3 ( L ) » ZKA3 ( L > ) - 

1  AKP3(1«»0«fY<L)»AK3(L)»ZKA3(L)))+BD3(L)*0.5*(AKP3(1.»0.»-Y(L+1) 

2  »AK3(L) »ZKA3<L) ) -AKP3 ( -1 . »0. r -Y < L ) r AK3 ( L ) » ZKA3 ( L ) ) ) )/EX3(L) 
TRAN3(L)  =  (AD3<L)  *AKP3  <**1«»0«»Y(L  +  1)»  AK3  (L)  »ZKA3(L)  )  + 

1  BD3<  L ) *AKP3 ( 1 . » 0« . -Y <L  +  1 ) » AK3(L > r ZKA3 ( L ) >  + 

2  ZKA3(L)*AR3(L)/(AK3(L)**2) )/(2.*EX3(L) ) 

GO  TO  50 

45  ZT3(L)=0,0 

TRAN3  <  L )  =  1 ♦ 0 
50  CONTINUE 

IF ( NSKAT « EQ ♦ 1 ) GO  TO  57 
IF (L«LE « 3 )  GO  TO  55 

ZK2 < L ) *ALOG (1«/<1«0-ZT2(L> ))/< (Y(L)-Y(L-l  )  ) ) 

ZK3 ( L ) -ALOG  <1./<1.0-ZT3<L)>)/((Y(L)-Y<L-1>) ) 

GO  TO  60 

55  ZK2  <  L ) =ALOG (1«/(1*0-ZT2<L) ))/<<Y(3)-Y(l))*.5) 

ZK3<L)=AL0G(1./<1«0-ZT3<L> ))/( (Y<3)-Y(l) >*.5) 

GO  TO  60 

57  ZK2(L )=ZKA2<L ) 

ZK3<  L ) =ZKA3  <  L ) 

60  CONTINUE 
M-L  + 1 

IF  <  L ♦ EQ . 1 )  GO  TO  65 

TAU2(L)=(ZK2(L)+ZK2<L-1) ) * ( YA ( L > - YA ( L-l ) ) *0 . 5+TAU2 ( L-l > 

TAU3<L')  =  (ZK3(L)  +  ZK3<L-1))*CYA<L)-YA<L-1 )  )  *0 . 5  +  TAU3  (  L-l ) 

65  EX2<M)=EX2<M-1)*<1.-(ZK2<H-1>*(Y<H)-Y(M-1))))+ZKA2(M-1)* 

1  EB2<H-1 >*< Y(H)-Y<M-1 > > 

EX3(H>*EX3<H-1 >*< 1 .-<ZK3(H-1 >*< Y<M)-Y<H-1 > ) ) ) +ZKA3 < M- 1 ) * 

1  EB3<H-1)*<Y<H)-Y<H-1)  > 

70  CONTINUE 
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U  fj  H 


c  SUM2=2.0*ZKA2< 1 )*EB2< 1 > +QS*TRAN2 < N > 

r  SUM3=2.0*ZKA3< 1>*EB3< 1 ) +GS*TRAN3 < N > 

GSW=GS*TRAN2<N)*TRAN3<N> 

DO  75  I=Nf2f-1 

C  SUM2=2.0*ZKA2< I >*EB2< I )*TRAN2< 1-1 )+SUM2 

C  SUM3=2.0*ZKA3< I >*EB3< I >*TRAN3< 1-1 >+SUM3 

GSW=GSW*TRAN2< I - 1 ) *TRAN3 < I - 1 ) 

75  CONTINUE 

CALL  RADWAL 

SUMR* < 1.0-REFT)*<SUMR-QSW)/3600« 

C  SUMR=< 1 .0-REFT)*<SUM2+SUM3-EB2< 1 >-EB3< 1 >  >/3600. 

RETURN 
END 


FUNCTION  AKP2 <SfPfXfAfB) 

AKP2=<1 .0+S*A/B)*< 1 «0+P*S*A/B)*EXP<2.0*A*X> 

RETURN 

END 


FUNCTION  AKP3<SfPfXfCfD) 

AKP3=<1.0+S*C/D)*<1.0+P*S*C/D)*EXP<2.0*C*X> 

RETURN 

END 


SUBROUTINE  RADCOMU) 

COMHON/BLKl/LSTfTSf AS f BS f TEMPA < 050 )  f RHOA < 050 ) f YA < 050 ) fEB2<050) f 
1  EB3(050) »TAU2(050) »TAU3<050) f ZK2 < 050 >  f ZK3 < 050 ) 

C0MM0N/BLK2/N fNP1fNP2fNP3f  NEQ  f  NPH  fKEXfKINfKASEfKRADfPEIfAMIfAMEf 
DPDXfXUfXDfXP  fXLfDXf INTO fCSALFAfHGSfC INF ( 10) fPREF  < 10) f 
PR (10) fP(10) fU(050)fF(10f050>  fR(050) , RHO ( 050 ) f ON ( 050 ) f Y ( 050 ) f 
DEN <050) * AMU <050) 

C0HH0N/BLK3/KRP.HACHfKRfINDI<10) fINDE<10) fSUHRfAKfALMGfBETAfTAUI 

1  F TAUE F VINF F DELTA F TW F TINF F QS F YL F UHAX » UNI NfFRf YIP F YEN fTRf HU F 

2  HF < 10) fGAMA< 10) fAJI < 10) fAJE< 10) fSC<050) fAU<050) fBU<050>  fCU<050) f 

3  A<10f050) fB<10f050) fC<10f050) fTEMPE<050) f TEMP < 050) fP0<050) f 
A  AMACH <050)fCP<050) 

C0MH0N/BLKA/NHfNN0fNNPfN0fNNfN02fNN2fNEf IND< 10) fSTOfAKSfRTfFT  f 
1  AMTfDRAGCfSTNOfTAUIUfQUfUGUfUGD 

DIMENSION  E12<200)fE13(200)fE12S<200)fE13S<200)fDIV<200) 
DIMENSION  AUX <  300 ) 

EXTERNAL  EBLT 
BEN=0 .075 
M  =  I--1 
DELT=0 • 1 
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1  S- TfMPA  (  M  )  -fOCl  T 
kHOS“DEN<H>*TR/TS 

7K2S=< 4. 37/0. 03281 )*< < RHOS/BEN ) ** 1 . 009 ) * < ( TS* 1 . 8E-04 ) **2 . 85 ) 
ZK3S=(4.985E-02/0. 03281 )*< < RHOS/BEN >** 1 .205) *<  <TS*1 ,8E-04>*» 
l  5.47) 

SUM2=0 . 0 
SUM3=0 . 0._ 

SUN2S=0 ♦ 0 
SUM3S=0 « 0 

CALL  EXPNI ( E2A2  »2»TAU2(M)»IER) 

CALL  EXPNI ( E2A3  » 2 »  T  AU3  ( M ) » IER  ) 

CALL  EXPNI < E2B2, 2» TAU2 <NP1 >-TAU2< M) » IER) 

CALL  EXPNI <E2B3 » 2» TAU3 <NP1 ) -TAU3 <H >» IER) 

CALL  EXPNI < E 12 <1 ) . 1 » ABS ( TAU2 < M ) -TAU2 ( 1 > )»IER> 

CALL  EXPNI < E13< 1 ) » 1 » ABS < TAU3 < M ) -TAU3 < 1 ) ) » IER ) 
TAU2S=(ZK2S+ZK2(M-1) >* (YA(H)-YA(M-l ) > *0 . 5+TAU2 < M-l ) 
TAU3S=<ZK3S+ZK3<H-1> )*<YA<H)-YA(M-1 ) > t0.5+TAU3(M-l > 

CALL  EXPNI (E2A2S»2rTAU2Sr IER) 

CALL  EXPNI < E2A3S 1 2 » TAU3S » IER) 

CALL  EXPNI <E2B2Sr2»TAU2(NPl)-TAU2Sf IER) 

CALL  EXPNI <E2B3S»2f TAU3<NP1 )-TAU3S» IER) 

LST=0 

EB2S=DML IN < EBLT » .05» ♦ 135» 1 » 300 r 0 . 0 » 0 . 001 » IER) 

EB3S=DMLIN<EBLT» .135f20. » 1 » 300» 0 .0 » 0 . 001 , IER ) 

CALL  EXPNI <E12S< 1 ) »1»ABS<TAU2S-TAU2<1) ) t IER) 

CALL  EXPNI <E13S< 1) » 1 » ABS < TAU3S-TAU3 < 1) ) r IER) 

DO  AO  L=2» NP1 

CALL  EXPNI <E12<L) » 1 » ABS < TAU2< H ) -TAU2 ( L ) ) » IER) 

CALL  EXPNI(E13<L)»1 » ABS ( TAU3 <  M ) -TAU3 (L))rIER) 

CALL  EXPNI (E12S(L) » It ABS < TAU2S-TAU2 ( L ) )»IER) 

CALL  EXPNI <E13S(L) » 1 »ABS<TAU3S-TAU3(L) ) rIER) 

SUM2=<EB2<L)*E12<L)+EB2<L-1)*E12(L-1) >*<TAU2<L)-TAU2<L-1) )*.5+ 
1  SUH2 

SUM3*<EB3<L)*E13<L)+EB3<L-1)*E13(L-1 ) ) * ( TAU3(L ) -TAU3( L-l ) )*.5  + 
1  SUM3 

SUM2S=<EB2<L)*E12S<L)+EB2<L-1)*E12S<L-1> )*(TAU2(L)-TAU2(L-1 ) ) 

1  * ♦ 5+SUM2S 

SUM3SMEB3<L)*E13S<L)+EB3<L-1 )*E13S(L-1 ) > * < TAU3 < L ) -TAU3 ( L- 1 ) ) 

1  ♦ . 5+SUM3S 

60  CONTINUE  - 

DIVRAD*2.*ZK2(M)*(E2A2*EB2(1)+E2B2*EB2(NP1)+SUH2-2.*EB2(M) ) 

1  +2.*ZK3<M)*<E2A3*EB3<1 )+E2B3*EB3(NPl ) +QS*E2B2*E2B3+ 

2  SUH3-2 . *EB3 < M ) ) 

DIVRAS=2.*ZK2S*<E2A2S*EB2< 1 > +E2B2S*EB2 < NP1 > +SUH2S-2 . *EB2S > 

1  +2,*ZK3S*<E2A3S*EB3<1 ) +E2B3S*EB3 < NP 1 ) +QS*E2B2S*E2B3S 

2  +SUN3S-2 . #EB3S ) 

ASM32.2*0.2162*DIVRAD/<RH0<I>*U<I))  ) 

BS=32.2*0.21A2*( <  DI VRAD-D I VRAS ) /BELT +  DIVR AD/TEMP < I ) )/<RHO<I >* 

1  U( I )*CP< I ) ) 

RETURN 

END 


i 
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SUBROUTINE  RADWAL 

COMMON/BLKI/LSTf TS»ASf BS»TEMF  A<050) »RH0A(050) »YA(050) »EB2(050)  r 
EB3<50)»TAU2<050)*TAU3<050)tZN2(050>rZK3(050> 

COMMON/BLK2/N»NP1# NP2 » NP3 » NEG » NPH » KEX » KI N » KASE » KRAD r PEI » AMI.AME, 
DPDX»XU»XD»XP»XL»DX»INTGfCSALFA»HGS»CINF(10>»  PREF  < 10 ) » 

PR  <10>?P(10)tU< 050 ) »F<10»050) »  R  <  050 ) » RHO ( 050 ) r QM  <  050 ) »Y( 050 ) » 

BEN ( 050 ) »  AMU  <  050 ) 

COMMON/BLK3/KRP r MACH » KR » INBI ( 10) f INDE< 10) » SUMR » AK » ALMG > BETA » TAUI 
rTAUErVINFf DELTA » TW » TINF f QS » YL » UMAX » UMIN f FR t YIP » YEM r TR » HU » 

HF(10) »GAMA< 10) »AJI<10>  »AJE< 10) » SC (050) »AU<050> » BU ( 050 > » CU ( 050 ) , 
A< 10 » 050) »B( 10 f 050) » C< 10 » 050 ) » TEMPE ( 050 >  t TEMP <050 >  » PO ( 050 >  » 

AMACH ( 050 ) fCP<050) 

DIMENSION  E2T 2(200) »  E2T3 ( 200 ) 

SUMR2=0 « 0 
SUMR3=0  «  0 
E2T2 ( 1 ) =1 . 0 
E2T3< 1 )=1 .0 

CALL  EXPNI (E3T2r3fTAU2(NPl)fIER) 

CALL  EXPNI <E3T3»3»TAU3(NP1)»IER) 

DO  15  K  =  2  r NP1 

CALL  EXPNI <E2T2< K) »2f TAU2(K) »IER) 

CALL  EXPNI (E2T3(K)r2»TAU3(K)rIER> 

SUMR2=(EB2(K)*E2T2(K)+EB2<K-1) #E2T2<  K-l )  ) * ( TAU2 (K) -  TAU2( K-l )  ) 

*  4 5+SUMR2 

SUMR3=<EB3(K)*E2T3<K)+EB3(K-1)*E2T3(K-1) >*(TAU3(K)-TAU3(K-1)  ) 

* • 5+SUMR3 
CONTINUE 

SUMR=2*EB2(1)-2.*EB2<NP1)*E3T2-2.*SUMR2 
+  2*EB3(1)-2.*EB3<NP1)*E3T3-2.*SUMR3 

REFT=0.0 

SUMR=  < 1 , 0-REFT ) *SUMR/3600 • 

RETURN 


SUBROUTINE  EXTRA 

C0MM0N/BLK2/N»NP1 ,NP2, NP3 » NEQ , NPH , KEX , KIN »KASE » KRAD » PEI , AMI ,AME, 

1  DPDX,XU,XD»XP,XL,DX, INTG,CSALFA,HGS,CINF< 10) ,PREF< 10) , 

2  PR  < 1 0 ) ,P< 10) »U<050> ,F<10,050> ,R<050>  »RH0<050> ,0M<050> ,Y<050> , 

3  DEN < 050 >, AMU <050) 

..  C0MM0N/BLK3/KRP»MACH,KR,INDI<10> »INDE<10>  » SUMR , AK , ALMG , BETA , TAUI 

1  »TAUE»VINF,DELTA»TUrTINF»QS»YL»UMAX»UMINrFR,YIP»YEM»TR»HW» 

2  HF(10) » GAMA < 10) ,AJI (10) , A JE ( 10 ) » SC ( 050 ) » AU<050> »BU<050) ,CU<050>  , 

3  AClOf 050) » B< 10 » 050) »C< 10,050) ,TEMPE<050> » TEMP < 050 ) , PO ( 050 > , 

4  ANACH(OSO) ,CP<050) 

COMM0N/BLK4/NH,NNO,NNP,N0,NN,N02,NN2,NE,IND<10) , STO » AKS » RT , FT , 

1  AMT,DRAGC,STNO,TAUIU,QU,UGU»UGD 
COMMON/K/KRI 

C0NM0N/6C/GC0N ( 050 ) ,AM0L(050) 

K I  N= 1 

NEQ-4 

KEX-2 

NPH-NEQ-1 

N*30 

NP1-N+1 

NP2-N+2 

NP3=N+3 

PINF=14 . 7 

TR=537 « 0 

UGC0N=1543 « 0 

AM0L1-29 » 0 

AM0L2-29.0 

DO  113  1*1, NP3 

302  AHOL ( I ) *1 • 0/ < <F(2» I )/AMOLl >+<F<3,I)/AhOL2) > 

113  CONTINUE 

DO  303  1*1 »NP3 

303  GCON  < I ) *UG CON /AHOL < I ) 

AMUl*1.2267334E-05 
AHU2*1 ♦ 2267334E-05 

DO  301  1*1, NP3 

301  AHU( I)*(F(2,I) *AMU1/AM0L1 )/(F(2,I) /AM0L1+0 ♦ 3535*  < ( 1 . +AM0L1/ 

1  AM0L2)**<-0.5) )*<  <l,  +  < AMU1/AMU2) **0.5* <  < AM0L2/AM0L1 > **0.25) )t* 

2  2. ) )+<F(3,I )*AMU2/AH0L2)/(F( 3,1) /AM0L2+0. 3535* ( (1 . +AM0L2/AM0L1 

3  >**(-0.5) )*<<1 .+<AMU2/AMU1) **0.5* ( < AH0L1/AH0L2) **0.25) )*t2.  )) 

DO  304  1*1, NP3 

304  DEN(I)-PINF*144./<GC0N< I >*TR> 

CP1*5990 • 0 

CP2-5990.0 
DO  306  1*1, NP3 

306  CP<I)-F<2,I) *CP1+F  <3,I)*CP2 

INTI AL  VALUE  OF  AUXILIARY  PARAMETERS 
TEMP  < 1 >  *TW 
DO  30Q  1-2, NP2 

TEMP<I)*<F<NH,I)-U(I)*U<I)/2»0)/CP<I) 

300  TEMPE< I >*TEMP( I ) 

CALL  DENSTY 


DO  351  1  ~  1 » NF‘3 

F'0(  I )  =DFN ( I )  *  <  TR/TEHF' <  I  )  )  *  ( GCON <  I  )  *TEMP  <  I  ) +0 . 5*U < I > /32 . 2 > /I 44 
AMACH<  I )  =U<  I  )/SQRT<  1 .  4*32 . 2*GC0N <  I  >*ABS<TEMP<  I  >  )  ) 

CALCULATION  OF  RADII 
CALL  RAD(XUfRd)fCSALFA) 

IF(CSALFA.EQ.O.OR.KRAD.EQ.O)  GO  TO  27 
DO  28  I=2»NP3 
R(I)=R(1)+Y(I)*CSALFA  . 

GO  TO  2? 

DO  30  1  =  2 » NP3 
R  < I ) =R  < 1 > 

CONTINUE 

CALCULATION  OF  OMEGA  VALUES 
OM  < 1 ) =0  « 

OM  <  2 ) =0 « 

DO  4?  I=3t NP3 

OM<I)=OM<I-1)+.5*<RHO<I)*U<I)*R( I>+RH0(I-1 > *U( 1-1 > *R ( 1-1 ) )* 

< Y< I )-Y< 1-1 )  ) 

PEI=ON  <  NP2 ) 

DO  5?  1  =  3  »  NP1 
0M( I )=0M( I )/PEI 
0M(NP2)=1 . 

OM ( NP3 )  =  1 . 

IF (NEQ«EQ* 1 ) RETURN 

DO  69  J*1,NPH 

IF(KEX«EQ«1)INDE(J)=1 

IF(KIN.EQ.i)INDI<J)=l 

CONTINUE 

RETURN 


